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ABSTRACT 


■  >  1  11  '  -y-'-k.  KuiK •*.  ^ 

This  dissertation  formulates  and  reports  upon  ihe  implementation 
of  a  numerical  system  for  the  solution  of  hydrodynamics  and  radiation 

diffusion  as  a  multi-material  problem  in  one  dimension.  A  parametric 

d-  I-  ‘  h,'  ■■ 

system  is  developed  in  which  the  program  parameters  may  be  dynamically 
altered  and  studied  as  to  their  worth  and  effectiveness.  The  system 
is  designed  specifically  for  use  within  an  interactive  man-machine 
environment  wherein  the  user  becomes  an  integral  part  of  the  final 
solution. 


vi 


1. 


d. 


INTRODUCTION 


The  largest  and  fastest  computers  have  always  been  used  for  the 
solution  of  partial  differential  equations,  especially  non-linear 
equations  which  are  used  to  describe  some  physical  phenomenon  in  time 
and  space.  Computer  programs  for  this  class  of  problems  are  large, 
and  their  creation  requires  a  joint  effort  of  many  individuals  over 
long  periods  of  time,  most  of  which  is  consumed  by  the  debugging 
process.  Once  such  a  program  has  been  developed,  its  use  requires  a 
great  deal  of  data  to  specify  the  desired  physical  system.  It  also 
requires  an  intimate  knowledge  of  the  workings  of  the  program  and  a 
vast  amount  of  intuition  and  experience  into  the  mechanics  of  the 
physical  processes  Involved.  Even  without  difficulties,  such  problems 
run  for  hours  at  a  time  on  the  most  modern  computers  in  the  typical 
batch  mode.  At  some  time  after  what  may  develop  into  days  and  weeks 
of  aborts,  restarts,  parameter  changes,  program  patches,  reconfigura¬ 
tions  and  the  like,  the  user  finally  acquires  several  edge  feet  of 
printed  output  and  perhaps  a  few  computer  generated  graphs  which  repre¬ 
sent  the  solution  to  his  problem.  He  must  then  examine,  plot  and 
otherwise  become  familiar  with  this  output  data  and  make  judgments  as 
t d  its  validity  and  applicability  within  the  constraints  of  the  system 
being  designed  or  simulated. 
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With  the  advent  of  multiprogramming,  time-sharing  and  real-time 
problem  solving  at  a  remote  console,  new  hardware  and  software  tools 
are  being  developed  to  allow  the  programmer  and  the  user  to  become  an 
active  part  of  the  checkout  and  running  process  of  a  program.  While 
this  document  reports  on  the  development  of  a  medium  size  program 
within  the  rudiments  of  such  an  interactive  environment,  its  primary 
emphasis  is  placed  upon  the  derivation  of  a  numerical  system  for  the 
solution  of  hydrodynamics  and  radiation  diffusion  as  a  multi-material 
problem  in  one  dimension.  The  system  is  developed  parametrically  in 
a  very  general  form.  Thus,  the  user  is  able  to  dynamically  configure 
the  system  into  a  form  best  suited  for  his  immediate  needs  through 
the  program  parameters.  The  technique  is  not  unlike  that  of  adjusting 
and  tuning  a  fine  piece  of  complex  mechanical  equipment.  It  also 
inherits  many  of  the  disadvantages  of  such  mechanical  systems,  primarily 
the  difficulty  of  dynamically  changing  the  program.  More  will  be  said 

on  this  and  the  requirements  for  man-machine  systems  for  these  types 

v 

of  computations  in  the  concluding  chapter. 

The  development  of  this  system  has  taken  place  over  the  period 
of  some  three  years.  During  this  time,  a  number  of  preliminary  computer 
programs  and  interactive  graphical  display  systems  have  been  written 
and  developed.  The  work  on  this  system  originated  at  Los  Alamos 
Scientific  Laboratories  in  New  Mexico.  It  then  moved  to  the  University 
ot  Utah  and  subsequently  to  Montana  State  University.  Computer  programs 
of  the  system  are  currently  operational  at  Los  Alamos  and  at  Montana 
State  University  through  remote  graphics  facilities  in  connection  with 
the  University  of  Utah. 
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Over  this  period  of  development ,  many  ideas  and  techniques  have 
been  explored,  tested,  accepted  and  rejected.  This  is  particularly 
true  with  respect  to  the  physics  and  the  subsequent  numerical  treat¬ 
ment  of  interface  conditions.  As  program  parameters  were  tested  and 
as  comparative  analyses  were  made  with  other  systems  and  solutions, 
both  analytical  and  empirical,  different  techniques  and  features  were 
incorporated.  These  changes  are  the  results  of  several  years  of 
experience  and  formal  education  not  only  in  the  areas  of  physics, 
mathematics  and  analysis,  but  also  in  computer  science. 

In  chapter  2,  the  difference  approximations  to  the  partial 
differential  equations  are  derived  and  the  complete  system  of  solution 
is  presented.  The  remaining  chapters  are  used  to  give  the  details  of 
the  auxiliary  calculations.  The  volume  and  mass  center  calculations 
are  given  in  chapter  3,  and  chapter  4  discusses,  the  time  step  selection 
procedure  and  the  associate  restrictions  and  control  parameters. 

Chapter  5  deals  with  the  material  properties  and  the  calculational 
aspects  of  the  various  thermodynamic  and  opacity  coefficients.  Chapter 
6  discusses  the  various  ways  of  specifying  and  calculating  the  source 
terms,  and  chapter  7  is  c  summary  of  calculational  results.  Chapter 
8  concludes  with  some  comments  on  future  research  areas,  particularly 
with  respect  to  the  man-machine  systems  alluded  to  above. 

This  chapter  is  concluded  with  definitions  of  the  symbols  and 
units  of  measure  used  throughout  the  remainder  of  the  text. 


1-1  Nomenclature 
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Symbol 

Description 

Unit  of  : 

A 

area 

2 

cm 

b' 

gas  constant 

j/Kv-j 

e 

internal  energy  per  unit  mass 

j/gm 

E 

internal  energy 

j 

m 

mass 

gm 

p 

pressure 

j/cc 

r 

radius 

cm 

s 

source  per  unit  mass 

j/gm 

S 

source 

j 

t 

time 

sh 

u 

velocity 

cm/ sec 

V 

specific  volume 

cc  /  gm 

V 

volume 

cc 

X 

space 

cm 

X 

•  • 

velocity 

cm/ sh 

X 

acceleration 

2 

cm/sh 

p 

density 

gm/ cc 

0 

temperature 

Kv 

<p 

fourth  power  of  temperature 

Kv4 

K 

opacity 

2/ 

cm  /gm 

I 

Rosseland  mean  of  the  mean  free  path 

cm 

V 


5 


1.2  Units  of  Measure 

gm  gram 

cm  centimeter 

j  jerk 

Kv  kilovolt 

sh  shake 

1.3  Physical  Constants 

a  radiation  density  constant 

c  velocity  of  light 

1.4  Conversion  Factors 

1  atmosphere  =10  ^j/cm  =10  kilobars 
1  j  =  1016  ergs 

7  o 

1  Kv  temperature  equivalent  =  1.16049  X  10  K 
-8 

1  sh  -  10  seconds 


.013732  j/cm-Kv4 
299.7925  cm/sh 


2. 


£ 


SYSTEM  OF  SOLUTION 
FOR 

HYDRODYNAMICS  AND  RADIATION  DIFFUSION 

A  set  of  difference  equations  and  a  system  for  their  solution 
Is  developed  for  hydrodynamics  and  radiation  diffusion.  The  momentum 
equation  is  differenced  in  a  natural  way  assuming  an  average  density 
at  the  interface  in  lieu  of  the  standard  area  over  mass  technique. 

Thus,  an  actual  pressure  gradient  is  computed  between  pressure  points 
calculated  at  centers  of  mass. 

From  the  beginning,  the  goal  was  to  difference  the  energy  equa¬ 
tion  in  terms  of  the  temperature  to  the  fourth  power.  This  approach 
was  selected  because  it  appeared  to  be  the  most-  natural  and  least 
complicated  in  contrast  with  the  more  traditional  differencing  schemes 
in  terms  of  the  temperature  or  change  in  temperature.  A  fully  parameter¬ 
ized  system  was  developed  in  very  general  terms.  This  permits 
detailed  studies  into  the  effects  of  time  differences,  interpolation 
and  extrapolation  functions,  smoothing  functions  and  the  like.  In 
addition,  an  iterative  procedure  is  employed  to  preserve  the  non- 
linearities  with  respect  to  the  energy  derivatives,  pressure  and  mean 
free  path. 

The  complete  set  of  difference  equations  together  with  boundary 
conditions  are  solved  in  a  well  defined  sequence  over  an  incremental 
unit  of  time.  The  solution  is  represented  by  temperatures,  pressures 
and  densities  as  functions  in  time  and  space. 
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2.1  Partial  Differential  liquations 

The  mass,  momentum,  and  ane  gy  conservation  equations  are  solved 
together  in  a  Lagrangia.i  system  where  the  motion  of  fixed  mass  points 
is  followed  in  plane  geometry.  The  substantial  derivative  forms  of 
these  equations  in  vector  notation  [lj  are: 


^  *  -P  (V'u) 


du 

p  ^  -  -?p. 


de  it,  „  *•  ,  ds 

p  dt  “  ~p(v'u;  *  v*q  +  o 


dt  ’ 


(2.1) 


(2.2) 


(2.3) 


where  the  viscous  pressure  and  gravitational  terms  have  not  been 
included.  The  Rosseland  radiation  diffusion  equation  [2], 

q  -  -  Y  •*> 

is  used  for  the  flux  term  and  the  standard  conduction  term  is  omitted. 

Equations  (2.1)  and  (2.3)  are  combined  to  give  a  system  of 
equations  for  hydrodynamics  and  radiation  diffusion. 

0  It  "  "7p*  (2.5) 


dv 

"ppdl 


lV<t, 


,  ds 
+  pdt 


(2.6) 
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2.2  Hydrodynamic  Difference  Equations 

Equation  (2.5)  1*  differenced  to  describe  the  motion  of  interfaces 
defining  the  constant  mass  zones.  Figure  1  illustrates  the  differencing 
scheme  where  the  center  of  mass  of  a  zone,  denoted  by  the  half  index 
i+*i,  givas  the  position  of  the  temperature,  pressure  and  density.  The 
Interfaces  delineate  the  mass  zones  and  are  referenced  through  the 
Integral  indices  1  and  1+1.  Note  that  these  Interfaces  are  fictitious 
interfaces  arbitrarily  constructed  to  obtain  difference  equations  repre¬ 
sentative  of  the  actual  differential  equation  given  by  equation  (2.5). 
Real  Interfaces  between  material  types  are  maintained  and  Included 
within  the  difference  scheme. 


Figure  1 
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-ccul.ration.  end  spatial  coordinate.  are  calculotjd 

espllcltly  at  each  intorf.ee  through  th,  folWag  ..,u.r.c.  of  e<!uatlo„s 
i/n  .  ,,n  n  n 


Vn  +  vn 
1  "l-n,  +  “>1_ia 


n  n 

Pl-%  ~  Pi. 


(2.7) 


^  .  -n-*s  +  ..n  Un^:  +  un~^  . 
i  i  £  2  ‘ 

xn+l  „  xn  +  ^  . 

whero  the  buPBr!Jcript  n  danotes  the  cln«  t.c"  wich 

tn+1  -  £n  + 


(2.8) 


(2.9) 


(2.10) 

Motion  (2.7)  t8  the  on.-di„n.iOK11  diff.r.nco  fora,  of  .gu.tlon  (2.5) 
-hioh  describe,  the  motion  of  the  interf.c,  at  ^  b.twe.n  the  .d)ac.„t 

*’ace  coordinates  >n+1  and  *n+1  r  ,  n+1 

1+1’  the  now  voJunn  vj^  is  calculated 

taking  ln:o  account  the  dimensions  of  tha  zone  1-^. 

Pr°8‘'an  1S  CaP'lbl“  01  “‘“^'‘n*  the  voiuno  olenents  In  on. 
o'  three  geometrl.,:  plan.,  cylindrle.l  .„d  .p„.tlc.a.  In  pl„„, 

geometry,  th.  r.dlu.  of  th.  cont.in.r,  often  referred  t0  ^ 

-y  be  constant  or  chon,.  in  .  ll„..cly  continuous  „  dl|cp„tlnuous 
fashion  ..  .hown  in  figer.  2.  Thls  flexibility  fa  p<r. 

"U‘  °CtU*'  PhySlCal  “  *  — *  -re  accurately  uith  rasp.  , 

to  velum.,  a„d  „r..s,  intruduclng  .  teo-di0.„,i„„.l 

■  ‘  0.,  that  th.  program  a.,u»cs  that  th.  u,«r  will  proceed  „lth 

caution  and  good  physical  intuition. 


Figure  2 


T. le  details  of  these  calculations  will  be  presented  later  in 
chapter  3. 

Given  the  new  volume,  an  avaro£«  density 


n+1  nl+4  1 

i+*S  ,n+l  n+1 

li*  vi* 

and  change  in  specific  volume 


£v 


n+*j 


vn+l  _  ,n 
1+*1  %L-*t 

®ul. 


(2.11) 


(2.12) 


are  calculated  for  use  within  the  diffusion  equation. 

The  motion  of  the  boundary  Interfaces  is  dependant  upon  boundary 
pressure  conditions  of  the  form 

P  *  P(0.  (2.13) 

where  p(t)  Is  given  by  a  user  supplied  subroutine.  Thus,  at  spntial 
position  i  ”  0, 


1  «n 


-  *S>/40 


(2.14) 
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and,  likewise  at  i  •*  N, 

n  ,  n. 

Vi  ‘  P<‘  >K 


••n  n 

XN  VN-*a  n  n  w , 

(XN  -  xK-V/6H 


(2.15) 


The  velocities  and  space  coordinates  are  calculated  using  equations 
(2.«)  and  (2.9).  The  change  In  volumes  at  the  boundaries,  Av”**1  and 

n+jj  0 

avn  atu  then  calculated  over  the  respective  intervals  (xj,  xn+1)  and 
r  n  n+1.  , 

'xjj  ’  XN  1  ort  ul  Co  obtain  tlic  work  done  by  the  boundary  pressures 
upon  the  adjacent,  zones.  These  energy  teres 


ASp 


n+l;  -  P(V0  <  ■  ,  and 


dSpf’  -  -p(cnj„  dV"41* 

i'  n 


(2.16) 

(2.17) 


are  simply  added  to  the  source  terras  of  the  firat  end  last  zone,  res¬ 
pectively. 

2.3  Radiation  Diffusion  Difference  Equations 

The  integration  of  equation  (2.6)  over  r.  homogeneous  volume  element 
and  the  application  of  the  divergence  theorem  yields 


dc  dv  ac  / — 

-  .  -ap  -  .  r  J  A, 


AVi  •  n  dA  +  m  ~ 
at 


In  ono  dimension 


Vi  •  n 
therefore, 


ex 


At  iL 1 

\  •  “7^  * 

>  p 

dv 

dV  * 

“J 

t  &£  /Til  HA  4 .41 

3  /A  )x  U  +  dT  • 


(?.18) 


(2.19) 


(2.20) 


with  tin  alg-  being  determined  by  the  direction  of  integration. 
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The  titno  advancement  schema  is  parameterized  bo  that  the  range 
from  oxplicit  through  fully  implicit  differencing  can  be  examined. 
Tills  generalized  differencing  permits  a  detailed  analysis  of  the 
numerical  processes  including  stability  questions  [3].  Also,  direct 
comparisons  between  the  explicit,  implicit  and  time-cer.tared  schemes 
can  be  mado  dyna-rlcally  with  the  same  program  on  the  same  physical 
problem. 

The  finite  difference  fora  of  equation  (2.20)  for  a  typical 
Interior  zone  is: 


n+1 


k  k+1 

Ui  ‘  Vs 


+  *i+i4  -  *144 


- “i(v> ' *.-j} 


n+S 

i+4 


where  the  superscript  k  denotes  tha  Iterative  value,  and 


k  r1 

— —  1  *  *u, 

M  jl*jj  1 


tl  «••«> 


n+1 


1+4 


k-1 


(l"Uj) 


(C 


(2.21) 


(2.22) 


(2.23) 
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V  'nt‘i 


Av 


i+>i 


Ai_ 

Av 


i  n+l 


Ac 

A* 


'  ,nf?  k/  \r.+l 

fcl  .  "  S  fe  ^>1  +  (l-u9) 


k-1 


fl+*S 


Ac 

Av 


ln+1 


i+*i 


k  n+*j  _  ,  ^  n+1  nl 

pi^-isp 


n+l  k 

Pi+»}  "  u3 


k-1 


n+l 


p  («,s)|^Ml-Uj)  P1+i,  , 


n+l 


aC  ,_*£ 

3  ^ 


n+l 


At 


n+*j 


(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 


iS14-1  *  pf  “s<«7.  (7.29) 

1  1+S 

Th«  function*  —  («.o),  ~  (6,o),  p(9,o)  and  th*  ratio  ar«  piven  in 
aectins  5.1. 

The  Bean  free  path  at  a  t/plcal  intarijr  intarfaca. 


*rl 


lT(, 

i-v 


tn+1 


*1+V  °i-y  °nV 


k-1 


+  (l-w4) 


rn+1 
i  ' 


(2.30) 


la  a  vcig  ted  average  rafla.tlva  of  cha  conditiona  near  tha  interface, 
the  detail*  of  which  are  given  ir.  section  5.2.  Th*  source  tens  consists 
of  a  space-dependent  part,  Pf,  often  referred  to  as  th*  power  factor, 
and  a  time-dependent  part,  AS(t).  Both  are  discussed  in  section  2.5. 
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Tlio  weights,  ,  u^(,  And  are  used  as  smooching  parameters 
to  dampen  out  tlio  discontinuities  introduced  by  linear  interpolation 
within  the  Equation  of  State  and  Opacity  tables,  A  typical  value  for 
these  weights  is  4,  Although  it  is  to  the  user's  advantage  to  use  a 
Vtt’uo  of  unity  it  the  tabular  values  are  smooth  or  if  a  continuous 
function  is  being  used  In  lieu  of  the  tables.  This  will  usually 
decrease  the  number  of  iterations  required  for  convergence.  On  the 
other  hand,  more  dampening  may  be  required  at  tinea  In  order  to  get 
convergence  at  all. 

The  finite  difference  equations  for  the  boundary  2ones  are  derived 
assuming  that  the  temperature  boundary  conditions  used  are  of  the  form 


u£  +  S  L  * 
ax 

with 


(2.31) 


|a|  +  |8  4  0. 


(2.32) 


In  terms  of  th«  derivative,  equation  2.31  becomes 


Dx 


'0 

0  -  *0  "  xq)^q 


*± 

bx 


‘  Vr.-S 


+  ««  <** 
»*  it 


xk-4)/6n 


(2.33) 

(2.34) 


at  each  of  the  boundaries 
first  or do i  approximation# 
equation  (2.31): 

{o  •  -  %  (\  •  V^o 


,  respectively,  oy  making  use  of  the  follcvlng 
for  «  at  the  boundaries  in  conjunction  with 


"nd  ♦»  ’  V«,  *  <*k  •  (2-3S) 


' 


.15 

If  6  -  then  it  is  not  exactly  corract  to  spoak  of  $  at  cho 
boundary  for  it*,  position  is  the  same  distance  from  the  boundary  as  la 
tho  adjacent  zono  temperature  but  in  the  opposite  direction.  The 
technique  is  often  referred  to  as  the  introduction  of  "fictitious  mesh 
points  about  the  boundary.  The  schurau  most  nearly  approximates  tho 
gradient  at  tho  boundary  as  if  there  were  an  additional  zone  extending 
out  from  the  bou\dary  whose  temperature  (flux)  and  pressure  profilos 
are  those  as  given  by  equations  (2.31)  and  (2.13).  If  6  -  l,  then  the 
boundary  condition  is  said  to  apply  exactly  at  tho  boundary. 

The  user  is  required  to  supply  a  subroutine  to  calculate  a,  fi, 
and  y  at  each  boundary.  They  are  either  constants  or  functions  of 
time  nnd/or  thermodynamic  variables  of  the  boundary  zones.  For  example, 
if  a  -  1,  fi  "  3,  and  y  "  t3(c)]\  then,  tho  boundary  condition  will  ba 
a  temperaturo  profile.  If  5  •  1,  then  this  profile  is  imposed  directly 
upon  the  boundary  interface.  This  typa  of  boundary  condition  impliaa 
that  the  b  undury  zonu  is  adjacent  to  a  raservoir  of  heat  that  can  supply 
or  absorb  energy  to  or  from  the  boundary  zone  respactively ,  yet  maintain 
a  predetermined  temperature  regardless  of  what  the  boundary  zone  or  any 
other  zone  in  the  problem  does.  This  boundary  condition  is  often  used 
to  facilitate  coupling  the  results  of  another  computer  program  to  this 
program. 

In  the  same  manner,  if  o  »  0,  3  -  1  and  y  ■  0,  then  no  energy  will 
bo  allowed  to  cross  the  boundary.  Thus,  a  perfect  insulator  at  the 
boundary  an  be  easily  specified,  or  a  syaaacrical  problem  can  ba  sol /ad 
~r  fently  by  doing  only  half  of  it.  In  othar  situations,  it  is 
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desirable  to  specify  the  flux  or  the  rate  of  energy  into  the  boundary 
zone.  This  can  be  done  by  setting  e  •  0,  6  - \  and  y  «*  F(t)  where 
F(t)  is  a  flux  profile.  The  mean  free  path,  A,  at  the  boundaries  is  a 
function  of  the  boundary  temperature  and/or  boundary  zone  temperature, 
and  the  boundary  zone  density.  Again,  the  details  are  deferred  to 
section  5.2. 

The  derivative  expressions  given  by  equations  (2.33)  and  (2.34)  are 
substituted  into  equation  (2.20)  to  obtain  a  finite  difference  approxima¬ 
tion  to  the  energy  equation  for  each  of  the  boundary  zones.  The  complete 
system  of  equations,  cne  for  each  zone  temperature,  forms  the  tri-diagonal 
system: 


„  n+l  .  „« 
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n+l  k  n+l  k+1  n+l 


-  a  u 
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(2.36) 
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k  n+1 
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(2.44) 


2.4  Implicit  Solution  for  Radiation  Diffusion 

The  solution  of  the  tri-diagonal  system  represented  by  equations 
(2.36),  (2.37),  and  (2.38)  is  an  adaptation  of  the  Gaussian  elimination 
scheme  for  diagonally  dominant  systems  [3].  Note  that  when  a=0,  that 
the  (j>'s  can  be  obtained  directly,  but  this  fact  is  not  exploited,  and 
the  following  scheme  is  a  generalization  for  all  values  of  a. 

Form  the  terms 
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2.5  Time  Advancement  Procedure 

Figure  3  is  a  simplified  flow  chart  which  illustrates  the  order 
in  which  the.  combined  hydrodynamic  and  radiation  diffusion  calculations 
are  done  during  one  time  step.  Note  that  the  hydrodynamic  calculi tions 
are  done  first  and  are  based  upon  the  pressures  from  the  previous  time 
step.  The  radiation  diffusion  calculations  are  then  done  using  the 
new  volumes  resulting  from  a  change  in  position  and  length  of  the  res¬ 
pective  mass  zones. 


Figure  3 
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The  rime  .top  ,t„  i.  defermln.d  by  .  ...  of  ujtr 

StabUlty  “*“*«.  «•  «£  which  1.  dependent  upo„  t„. 

accciototioc.  Ch-pter  «  is  devoted  „mplet.ly  t#  .  #f 

these  time  step  restrictions  and  the  procedure  ue.d  to  determine 

0Ptl“““  tl,a  St  increment  for  the  current  stability  require¬ 
ments  end  the  deeree  of  resolution  end  accuracy  desired. 

Subsequent  to  selecting  the  time  step  .ire  end  calculating  the 

new  spatial  dimensions,  an  extrapolation  scheme  1.  introduced  to 

provide  initial  temperature  estimates.  Its  purpose  is  to  reduce 

the  number  of  iterations  through  a  linear  oraWn 

sn  a  linear  projection  of  the  previous 

temperature  change: 


°i-+^ 


(2.53) 


The  extrapolation  parameter  v  will  tvDir»n«  . 

v  win  t/picelly  vary  between  0  and  1  with 

0,  h  and  1  being  the  more  popular  values  v.,^u 

F  puxar  values.  Further  comments  with  exampl 

arS  8iV6n  ln  SSCtlon  7 -1  w*th  respect  to  its  use. 

Given  an  initial  temperature  estimate,  the  quantities 
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"*  then  '«•  «<*  *°n.  In  eh.  normal  „r  ulthc  jt  any 

smoothing  parameters.  Likewise,  at  the  lnt,rfaCf 


ices 


0  n+1 
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n+1 

V  n4)l  . 


0  n+1  f  n+i 

1  I  (Vy  l+ij*  pi-i3*  “i+^j  for  i-1 . N-l,  nnd  (2>58) 
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(2.59) 


All  Initial  values,  .sc.pt  eh.  sour.,  cerru,  ar.  now  .v.n.bl.  t. 
.ease  eh.  le.raelv.  process  fcr  r.dl.clon  Jlfl-u,lon,  ^  calculatlon 

for  eh.  source  terns  1.  complicated  sine.  eh.r.  .r.  four  way.  of 

introducing  energy  Into  the  problem.  Two  have  b..„  discussed,  th. 

imposition  of  boundary  temperature  (flus)  and  pr.s.ur.  profiles.  N„t. 

that  the  boundary  temperature  or  flux  profil.s  do  not  d.r.etly  r.sule 

in  a  sourc.  term  calculation.  In  this  c...  the  energy  transfer  1. 

implicitly  included  in  the  energy  agnation.  heaver,  the  net  flow  of 

energy  across  th.  boundaries  due  to  th.  imposition  of  temperetur.  .„d/0r 
flux  profiles  Is  simply; 
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Thfi  °Cher  CV0  “ource  cern>  forms  are  discussed  ir, 
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chapter  6. 
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Tho  procuss  ituratos  as  shown  in  figure  3  until  the  temperatures 
converge  to  a  stable  valuu  Hint  is  until  for  each  rone 


k+1  n+1  k  n+l 
U14lj  "  '  14'.. 

k-n+1 


£  t  , 


(2.62) 


given  some  u*uf  specified  epsilon.  Since  convergence  in  not  guaranteed 
there  is  an  uppur  Unit  upon  the  nuabrr  of  so  called  temperature  itera¬ 
tions.  Twenty  is  a  nominal  iiaximum  number  of  iterations,  but  rhe  number 
may  be  altered  by  the  user  to  fit  Ills  problem,  in  particular,  the  limit 
may  be  set  to  unity,  rjducing  the  calculations  to  u  non-iterative  pro¬ 
cedure.  The  user  may  also  rei uce  the  time  step  size  if  the  number  of 
iteration*  per  time  step  becomes  excessive  or  it  the  temperatures  do 
rot  converge  below  the  specified  epsilon  of  convergence.  In  addition 
ti  the  time  step  size  restriction  parameters,  several  of  the  other 
averaging,  smoothing  and  extrapolating  parnmeters  discussed  above  affect 
the  rata  of  convergence  and  \he  number  of  iterations  rnquired  for  conver¬ 
gence.  The  ability  to  change  thesa  parameter!  dynamically  naka  it  pos¬ 
sible  to  give  the  user  the  answers  he  desires  in  a  most  optimum  way. 

It  vary  often  provides  the  only  •  eans  of  getting  past  a  particularly 
difficult  part  of  the  problem.  It  al«o  provides  the  means  of  getting 
over  particularly  stable  portions  without  excessive  computer  time.  Thus, 
the  program  need  not  run  always  with  these  parameters  set  at  worst  caso 


■\cvel6 . 
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Thd  1* «t  atop  in  tha  tJ*e  advancement  scheme  la  cha  updating  of 
tha  anargy  ter-s.  These  tjrms  ara  calculatad  and  displayad  aa  a  moana 
of  i.aeplng  track  of  anargy  conservation.  Thus,  tha  following  total  and 
par*  il  sums  are  maintained: 


n+S 
Ji-Hj 
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(2.65) 


For  anargy  conservation,  it  should  be  tha  caaa  at  all  tlsiaa  that 

“n+1  *  ^+l  *  ^  •  (2.66, 
This  will  never  ba  tha  case,  however.  Just  due  to  normal  truncation 

«nd  roundoff  .rrora  within  the  computer,  thaaa  a urn.  are  prad.atined  to 
not  add  up  correctly.  It  la  tha  relative  diff.rence  which  la  import¬ 
ant  for  tha  user  to  monitor.  Ha  can  affect  control  over  thia  diffaranca 
primarily  through  the  convergence  parameter  and  tha  time  a tap  raatrictioi 
factora  as  givan  in  chatter  4.  The  energy  coneervation  and  computer 
time  are  o  rnal  trade-offs.  lr  haa  been  my  experience  thot  tight 
energy  conservation  will  not  always  give  significantly  different  results. 
It  seems  to  ba  Important  only  at  certain  times  within  tha  problem.  Many 
of  chase  are  known  a_£riori  with  respect  to,  aay,  the  characteristics  of 
the  users  energy-time  profile  imposed  upon  a  boundary  or  a  set  of  aonec. 
In  the  main  however,  the  user  will  only  tighten-up  certain  parameters  as 


the  need  arises. 


3. 


ur 


ARU,  VOLUME  AMD  CENTER  OF  MASS  CALCULATIONS 

Attar  th*  Bccalnration.  velocity  and  space  terns  arc  advanced  In 
tine,  new  zone  volume,  ere  calculated  in  order  to  determine  tit*  changes 
in  epecll.c  volume  over  the  « tme  tine  interval.  These  changes  are 
indicative  of  the  amount  of  energy  converted  into  kinetic  energy  due 
to  preesurt  gradient*.  Indeed,  it  ia  of  great  importance  to  monitor 
the  aawunt  ot  energy  converted  and  its  rate  of  conversion  into  or  frost 
kinetic  energy.  This  energy  appears  as  shocks  or  shock  waves  and  their 
characteristic  pressure  peaks  can  be  dynamically  displayed  upon  a  CRT 
terminal.  The  user  is  then  able  to  monitor  the  motion,  distribution 
and  effects  of  such  shocks  In  direct  relation  with  concurrent  displays 
of  temperature  and  density  terms.  Of  special  interest  are  the  changes 
in  volume  at  the  problem  boundaries.  It  is  here  that  external  forces, 
in  the  form  of  boundary  pressures,  pump  kinetic  energy  into  or  out  of 
the  problem.  ln  addition,  the  motion  of  zone  boundaries  requires  that 
new  zone  mass  centers  turn  be  calculated  from  which  finite  difference 
terms  can  be  forced  to  describe  the  traneport  of  energy  from  one  zone 
to  the  next  by  the  diffusion  of  redietion. 

ihc  purpose  of  this  chapter  Is  to  give  the  computational  details 
usad  for  aeterining  the  volume  end  mass  centers  for  each  zone.  Each 
of  the  three  geometries  is  discussed, 

3- 1  Plano  Geometry 

In  piano  geometry,  the  space  dependent  radius  adds  to  the  coe- 
P*i  it  of  the  ares,  volume  and  cantor  of  meas  determinations.  First, 
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l«t  m  point  out  that  tho  center  of  hoss  w«  wish  to  determine  It  not 
^ln  kcnorni)  equivalent  to  tht  cantor  of  gravity.  Instead,  it  is  the 
point  on  tho  ax.U  which  divide*  the  r<  ne  in  half  with  ratpict  to  It* 
voluao  and  tonuequently  its  m»*  tinea  it  la  attuned  to  have  a  uniform 
density.  Secondly,  note  that  equations  (2.7)  and  (5.36)  require  such  a 
division. 

Unfortunately ,  else  hat  not  allowed  an  exact  division  for  the 
cate  In  which  the  radii  changes  linearly  with  respect  to  tpace.  Instead, 
the  center  of  gravity  was  calculated  as  an  approxiaation  for  the  center 
of  these  sections.  I  refer  to  such  sections  as  truncated  frustruas  or 

trapezoids  of  revolution. 


The  user  must  provide  a  description  of  the  contnlner.  This  Is 
dona  Independently  of  the  initial  spatial  description  of  tho  aatarlal 
within  In  the  fora  of  a  space  profile  of  radii  and  space  pairs 

\  ■  V' J '  *• 2 . ai) 


such  that 


*n  £  xn 

PJ  PM 


(3.2) 


and  r  >  0.  / 

Pj  < 

Tha  volut:a  of  a  zone  delialted  ay  Interfaces  at  x”  and  x"  at  soae 
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In  the  special  (and  cost  often  encountered)  case  for  which 
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If  the  restrictions  given  above  in  aquations  (3.12)  or  (3.13)  cannot 
ba  satisfied  because  pari  ov  all  of  tha  zone  lias  outaida  of  the  pipe 
profile,  then  the  formula  above  is  used  squally  wall  to  extrapolate  for 
tha  interface  radii.  In  these  cases,  either  equation  (3.4)  or  (3.14) 
apply  depending  upon  whether  more  than  one  (the  last)  pipe  section  is 
involved.  Care  must  be  token  to  insure  that  r“  >  0  for  all  i  -  1,  ....  N. 

The  zone  center  is  then  given  by 
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if  restriction  (3.13)  applies. 

The  area  of  e  typical  Interface  is  siaply 
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3.2  Cylindrical  and  Spherical  Geometry 

Hie  calculations  for  cylindrical  and  spherical  geometry  are  much 
simpler  than  that  for  plane  geometry  to:'  there  is  no  analog  to  the  flexible 
cross  section  allowed  therein, 
for  cylindrical  geometry 

A A  -  2*Xj  h  # 


(3.23) 


V  ■  f 

l-Ha 


and 


<*m>2  -  <^)2]  h . 


,  n  .2 


♦  <x”>2  i/ 


„ 

x1+ll-  A 

where  h  Is  Che  height  of  the  cylinder  as  supplied  by  the  user. 
For  spherical  geometry 


,n  .  ,  n.2 

Aj  •  4i?(Xj) 
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(3.24) 


(3.25) 


(3.26) 


i+S 


...  !.n  .3  ,  n  *  3 

4/3"  |(x1+1^  "  (-V 


(3.27) 


and 


Ci-^ 


-  s  ^  i\ 


>3  + 


,  n\3 

(x^J 


1/3 


(3.28) 


3.3  Change  in  Volume  at  the  Boundaries 

As  required  in  equations  (2.16)  and  (2.17),  the  change  in  volume 
at  each  boundary  is  calculated  In  a  manner  not  unllko  that  glvan  above. 
Since  the  boundary  volume  element  is  regarded  an  a  change  in  volume, 
the  sense  or  sign  of  the  change  is  important.  Note  that  the  volume  at 
boundary  i  ■  0  must  be  calculated  over  the  interval  from  Xq  to  Xq  • 

If  Xq+1  <  x£,  then  the  interval  is  thought  of  as  having  a  nagative 

niL 

length  if  order  to  determine  the  correct  sign  of  AVq  •  The  same  care 

and  consideration  is  given  to  the  other  boundary.  Note  that  if 

AVn  *  .  0  and/or  aVn+?s  <  0,  then  work  is  done  on  the  system. 

0  n 


TIME  STEP  SIZE  DETERMINATION 


A  number  of  restrictions  are  placed  upon  the  time  step  size  for 
purposes  of  stability,  resolution  and  convenience.  Each  is  discussed 
in  detail  with  respect  to  their  effect  upon  the  solution  and  each  is 
controlled  directly  or  indirectly  through  user  specified  parameters. 


4.1  Courant-Friedrichs-Levy  Hydrodynamic  Stability 
Criterion  [3j: 


4tn+li  *  At 


rA.xn+1 
n+^i  min  ‘  i+^ 

«a  y - — 


i 


lrfl 


(4.1) 


’i+.y 


This  restriction  prohibits  the  shock  front  moving  at  a  velocity  u 


from  passing  completely  through  any  one  zone  during  the  time  step.  Note 


that  and  u^+1  ate  not  known  at  the  beginning  of  a  time  step.  The 


n+^ 


zone  width  is  an  explicit  function  of  At  ^  and  the  shock  velocity  is  a 
function  of  the  resultant  pressure  and  density,  the  details  of  which  are 
given  in  section  5.1,  equation  (5.25). 

The  procedure  used  to  determine  the  Courant  time  step,  therefore, 


is  to  estimate  an  initial  value 


min 


ttf*  ■  c.  . 

c  fac  i 


(4.2) 


with  0  <  C  £  1,  and  then  check  after  the  initial  extrapolation  and 


each  subsequent  iteration  that 

r.  n+i 
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If  at  any  time  during  the  iterative  solution  it  is  found  that 

k-A,.n+^  „  .  .n-f^  . 

Atc  <  At  »  then  an  err°r  message  is  displayed  and  the  program  is 
put  into  a  waiting  mode.  The  user  is  then  expected  to  initiate  a  cor¬ 
rective  procedure  which  will  enable  the  program  to  continue  in  a  stable 
condition. 

The  corrective  procedure  may  be  nothing  more  than  reducing  the  time 

step  size  directly  or  indirectly  through  the  parameter  C.  and  then 

fac. 

backing-up  and  restarting  the  current  time  step.  It  can,  however,  involve 
a  detailed  examination  of  the  conditions  of  the  offending  zone  or  zones  in 
order  to  understand  the  mechanisms  creating  the  instability.  More  precise 
solutions  may  then  be  directed  at  the  particularly  sensitive  variables 
involved. 


4.2  Hydrodynamic  Zone  Increment  Change  Restriction 
n+1 


H  j.  max 
fac  i 


Ax 


1- 


Ax 


l+h 

n 

i+-a 


(4.4) 


This  is  a  restriction  on  the  fractional  change  in  the  length  of  a 
zone,  with  0  <  Hfac  <  1.  Its  primary,  purpose  is  to  minimize  the  discon¬ 
tinuity  in  the  distance  between  zones  and  the  volumes  of  each  zone  from 
one  time  step  to  the  next.  It  also  prevents  a  complete  collapse  of  a 
zone,  and,  in  particular,  the  inversion  or  crossing  of  adjacent  inter¬ 
faces.  In  addition,  this  restriction  eliminates  most  of  the  problem  of 
mismatched  zones  with  respect  to  their  mass  ratios.  The  problem  arises 
during  a  given  time  step  when  a  massive  zone  crushes  and  collapses  a 
much  less  massive  zone  before  the  latter  is  able  to  build  up  a  resisting 


33 


pressure.  Another  difficulty  with  adjacent  zones  of  widely  differing 
masses  occurs  when  an  average  opacity  is  calculated1  at  the  interface 
between  the  zones.  This  problem  is  discussed  in  detail  in  chapter  3. 
A  time  step  increment  satisfying  the-  restriction  given  by  (4.4) 


can  be  explicitly  calculated  once 


since 


.  n+1  n+1  n+1 

Axnh  •  xl*l  -  xi 


the  acceleration  terms  are  known 


n  -n 

Xi+1  *i 


+  -  i^)  At”* 

+  <*”  -  A  At”^  (Atn^  +  At"^) 


vi+l  ■  V  - 2 - '  (4-5) 

using  equations  (2.8)  and  (2.9).  Thus,,  a  maximum  value  of  tn+^  needs 


to  be  found  suen  that 
I a^(Acn+^)  2  +  b.  At^l  <;  Cl  , 

where 


(4.6) 


i  ,-n  -n 
^(x.  -  x 


i+1 


), 


a  At 


■  n-% 


+  (S^ 


.n-^  , 

Xi+1 


U  A  n 

H  Ax. 

fac  jl-Hs 


(4.7) 

(4.8) 

(4.9) 


for  i  -  0,  1,  ....  N-l.  a  solution,  At”**5,  of  the  system  represented  by 
(4.6)  is  in  general  a  quadratic  root  depending  upon  the  values  of  the 
coefficients  a.,  b.  and  c±.  Since  C±  >-0,  it  is  possible  to  show  that 
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viv 


-b,+(b^+4a,c,) 
i  i  i  i/ 


At 


n+%  min 

S3 

h  i 


2a . 

i 


~bi"(br4aiV'* 

2a. 


if  a^b^O 


if 

ai“0» 

b^O 

f 

a±>0. 

b^O 

if  < 

a±>0. 

bi<°, 

bi“4aic 

:i<° 

,v°- 

b1>0, 

bi+4aic 

/ai>0, 

bi<0” 

bi"4aic 

i*0 

it  < 

a.<0, 

i 

b^O, 

bi+4aic 

i<0 

^<0, 

b^O 

(4.10) 


The  user-  can  dynamically  change  the  importance  of  this  restriction 
through  the  parameter  H£m  ,  rezonlng  If  necessary  when  the  restriction 
is  too  greatly  localized  at  one  zone. 


4.3  Temperature  Change  Restriction 


At  £ 


At 


n+^s 
9  ' 


r  min 
fac  i 


9 


in 


i+% 


+  0 


I0i44s  0i4V 


At 


n-^ 


(4.11) 


This  restriction,  not  unlike  the  others,  is  added  to  reduce  the 
truncation  error  and  give  better  resolution  and  accuracy  when  desired. 
Usually,  0  <  Tfac  s  1,  and  the  parameter-  6p  is  a  reference  or  base 
temperature  above  which  the  restriction  is  to  be  applied.  Both  parameters 
may  be  set  dynamically  by  the  user  co  best  suit  his  requirements  at  any 
time  throughout  the  solution  of  his  problem. 


Specifically,  an  economic  trade-off  with  computer  time  consumption 
per  unit  of  problem  time  can  be  made  at  those  points  when  the  temperature 
and  pressure  gradients  are  such  that  a  first  order  approximation  remains 
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valid  over  a  much  larger  time  interval.  Note  that  the  three  time  step 
restrictions  discussed  up  to  this  point  will  automatically  increase  or 
decrease  the  time  step  increment  depending  upon  current  conditions 
internal  to  the  problem.  The  user,  however,  ha9  control  over  the  degree 
to  which  these  restrictions  will  be  applied. 

If 


At 


n-‘s 


A6 


i+Jj 


t  a*  0  , 


(4.12) 


meaning  that  the  temperature  change  occurs  irrespective  of  the  specific 
heat  of  the  zone,  then  the  temperature  change  restriction  on  that  zone 
is  ignored. 


4.4  Doubling  Restriction 


atn+ls  s  itf1-  -  2  « 

d 


n-Hs  „  . .  r. 


(4.13) 


This  arbitrary  restriction  prevents  the  time  step  size  from  becoming 
too  large  too  fast.  If  the  next  time  step  is  quite  large  with  respect 
to  the  previous  one,  corresponding  changes  in  .the  derivatives,  pressure, 
opacity,  etc.,  may  be  too  large  to  calculate  accurately  even  in  light  of 
the  iterative  scheme.  A  rapidly,  varying  time  step  from  one  time  step  to 
another  can  adversely  effect  the  hydrodynamics  through  equation  (2.8). 
Note  that  there  is  no  restriction  upon  how  small  a  time- step  can  become 
even  with  respect  to  the  previous  time  step  size. 
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4.5  Print  Restriction 


£  At 


n+^ 

P 


(4.14) 


At  user  specified  intervals,  the  program  prints  out  tha  currant 
space  coordinates,  temperatures,  pressures,  etc. ,  so  that  a  permanent 
record  of  the  problem  solution  may  be  retained  for  further  study  at  a 
later  time.  The  program  also  makes  a  hard  copy  of  the  image  on  the 
display  and  writes  out  onto  a  user  file  the  current  state  of  the  program 
so  that  it  may  be  restarted  at  that  time  if  so  desired  at  some  later 
time.  This  dump  feature  allows  the  user  to  run  a  problem  at  several 
sittings,  restarting  at  or  near  the  point  in  time  at  which  he  stopped 
previously.  It  also  facilitates  his  going  back  in  time  and  trying  a 
different  solution  path  by  changing  one  or  more  program  parameters  or 
variables. 


(4.17) 


if  th«ro  It  a  next  value  t  .  Also,  if 

tprtj+l  *  tprtj  +  mAtprtj  ’  (^‘18) 

t5>*n  Cprt  i8  sat  t0  CPrtj+i  *  80(1  th*  lnd,x  J  is  »dvanced  by  on*. 

Ih*  print  profll*  1*  used  to  obtain  print-outs  ov*r  certain  intervale 
of  Interest  at  "nice"  values  of  time,  usually  factors  of  two  and  five. 

Both  the  profile  and  the  next  echeduled  print  tine  may  be  altered  dynaa- 
ically. 

Since  the  print  times  will  normally  Interrupt  the  running  sequence 
with  a  short  time  step,  the  program  attempts  to  restore  the  time  step 
size  for  the  next  time  advancement  to  the  level  at  which  it  had  been 
running.  Specifically,  the  temperature  change  and  doubling  restrictions 

h«il  i 

will  use  At  ana  guide  in  lieu  of  Atn~*  .  The  other  reetriction 
criterion  are  used  ae  stated  after  a  print  cycle. 

A. 6  Maximum  Restriction 

Atn+4s  S  min(Atm  .  Atfl)  .  (4.19) 

Once  in  a  while  it  le  convenient  to  set  a  limit  on  the  size  of  the 
time  step.  It  is  particularly  helpful  when  the  user  is  attempting  to 
follow  a  phenomenon  which  le  not  being  controlled  automatically  by  any 
of  the  other  time  step  limiting  procedures. 

A  one  time  only  maximum  may  be  set  through  At#.  It  is  automatically 
reset  to  infinity  for  the  next  time  step.  There  is  an  automatic  back-up 
scheme  which  may  occur  if  the  temperatures  don’t  converge.  When  this 
back-up  occurs,  Atg  is  set  to  one  half  the  value  of  Atn44*  and  the  calcula¬ 
tions  for  that  time  step  are  restarted.  An  overall  maximum  Atm  may  be  set 
and  remains  as  set  until  reset. 
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Not*  that  whil*  At  ^  may  b*  altered  directly  by  user,  It  ia  not 
effective  unleae  he  also  schedules  e  back-up  since  At04*1  is  recomputed 
at  the  beginning  of  the  next  time  step.  The  back-up  procedure  does  not 
restart  the  time  step  calculations  at  the  beginning  of  the  advancement 
procedure  as  illustrated  in  Figure  3,  but  instead,  restarts  with  the 
selection  of  the  time  step  size.  The  accelerations  and  previously 
discussed  time  step  size  restrictions  remain  valid  during  a  back-up. 


4.7  Hold  Time  Restriction 


4tn+1  s  atn<Ja  -  t  -  t" 

W  W 


(4.20) 


It  is  often  useful  to  set  a  time  t  at  which  the  user  wishes  to  put 
the  program  in  a  hold  or  wait  condition.  This  feature  insures  that  he 
will  be  able  to  put  into  effect  changes  at  certain  specified  problem 
times  crucial  to  the  overall  problem  solution. 

If  tv  <  tn,  then  the  restriction  does  not  hold  and  the  program 
automatically  sets  Aty  to  the  machine  relative  Infinity. 


5. 

MATERIAL  PROPERTIES 

Equation  of  state  and  opacity  data  muet.  La  provided  for  each  speci¬ 
fied  material.  These  data  are  available  in  tabular  fom  on  a  mass  storage 
device  and  are  automatically  read  into  tha  program  whan  called  for  by 

noma  during  tha  input  phase.  The  user  may  also  include  his  own  special 

purpose  tables  or  routines.  Sometlmas  tha  notarial  properties  con  be 

calculated  from  a  set  of  parametarlcad  equations.  This  may  ba  tha  case 

when  a  problem  is  run  for  which  there  is  an  analytical  solution  to  chock 
with. 

This  chapter  will  limit  its  discussion  to  the  standard  materiel 
property  tablss  and  ths  procedures  which  are  used  to  calculate  tha 
required  thermodynamic  quantities. 


5.1  Equation  of  State 


Tha  gasaous  aquation  of  atata  tables  era  organised  by  density  end 
temperature  within  each  material  type.  For  each  tabular  density  in  tha 
form  lnvj,  there  la  a  sequence  of  tempareturi,  internal  energy  and  pro¬ 
portionality  triplets,  (8,  *  ,  b').  ...  ,  such  that: 

•  J 

p£  “  V  *  (5.1) 


lnv.  <  lnv,  <  ...  <  lnv.  , 


(5.2) 


9J.l  '  ®j,2  *  *•*  '  ®j,k  * 


(5.3) 


with  J,  K  :  l.  Note  that  there  is  no  requirement  that  the  teblss  be  the 

same  sl2e  from  one  material  type  to  another,  or  that  8.  .  •  8. .  for 

j ,k  J+l,k 

any  j,  k  within  the  table.  A  typical  maximum  table  sire  would  be  J  •  20 
and  K  -  30. 
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The  material  internal  energy,  t*,  include*  the  kinetic  energy  of  the 
free  pertlcle*;  the  dissociation,  ionisation,- and  excitation  energy;  but 
not  the  redletlon  energy  given  by 

cr  -  *4v.  <5‘*> 

The  equation  of  etete  teblee  ere  genereted  by  other  computer  program* 

[4],  and  there  are  many  table*  in  existence  for  the  common  element*, 
compounds  and  mixture*. 

Glvsn  some  temperature  and  density,  (6»p)»  the  energy  derivetive* 
and  proportionality  constant  are  approximated  from  these  tables  through 
the  following  interpolation  formulae: 


lnv  - 

if  (e’o)  •  isr- 


lnvi  [Cgi+1.*+1  gin,t  _  Jj Jsti - 5iiil 

“  lnVj  \0J+1.4+1  '  ®J»*  ®J.k+l  "  J 


J+l 

+  gt,k^1  8^*-“  +  4a63v, 
9J,k+l  "  ®j.k 


(5.5) 


~  (e.o)  -  pe2 


lnv  -  lnv^ 

lbl+i.t+i -  bj+;,t . . 

.  b.i,k+i  ~  bi.-k 

lnvJ+1  -  lnVj 

®J+1,  4+1  "  ®J+i.i 

®j,k+i "  ®j,kJ 

®J,k+l  '  ®J,k 


(5.6) 


lnv  -  lnv. 

"loVl  -'i%  (bj+l  ~  bP  *"j  • 


(5.7) 


where 


bj  -  bj.k  ♦  <*  -  9-  -> 


bi.lc»l  ~  bl.k 


(j.« 


*1 


bj+l  “  bj+l,l  *  <0 


o  i  '  N+i.t 

J+1>‘  6j+i,m  - 


assuming  thsc 


InVj  £  lnv  £  lnVj+1  with  1  £  J  <  J, 
®J,k*  0s:  6j,k+l  withl£k<K, 


J+1.4 


£  0  £  8 


J+l,t+l 


with  1  £  £  <  X. 


(5.9) 

(5.10) 

(5.11) 

(5.12) 


If  ons  or  raors  of  tha  abovs  conditions  given  by  (5.10),  (5.11)  and 
(5.12)  do  not  hold,  than  equations  (5.5),  (5.6),  (5.7),  (5.8)  and  (5.9) 
are  used  as  extrapolation  formulae.  The  evils  of  this  extrapolation 
scheme  have  shown  up  often  upon  the  graphics  scope  in  the  fore  of  non¬ 
convergence  and  dlscontiuulty  ir.  successlva  values.  It  is  hoped  that 
through  the  graphics  monitor,  •  new  extrapolation  scheme  or  limits  upon 
the  existing  schema  can  be  established  which  will  aid  the  user  in  bring¬ 
ing  his  problem  up  to  those  temper etures  and  dansitles  which  are  within 
the  bounds  of  the  tables.  The  user  must  often  coaproalse  with  storage 
limitations  and  put  in  tables  which  cover  only  a  particular  area  of 
interest  with  respect  to  teaperatura  and  density. 

Equation  (5.5)  gives  an  approximation  for  ,  the  apacific  heat. 
Ihe  tables  are  constructed  for  linear  intarpolation  on  8;  therefore, 

equation  (5.5)  in  terns  of  $  is  unsatisfactory  for  an  approximation  to 

it 

■j-T  .  An  acceptable  scheme  makes  use  of  the  relation 

d  C  ^  cj_£  _d_9 

H  "  30  9+  * 


(5.13) 
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with 

li  M 

3#  *  l*  * 

Therefore, 

k/i£,»  >|n+^  k/M,«  .l"4'1  k  U >' 

44  JlW  I40  "’ll*  U  144, 

where 


(5.14) 


(5.15) 


k  foil  ■* .  !SsLj£& 


l-C +  *W  ([W 4  M) 


(5.’«) 


Note  th.it  this  approximation  avoids  problems  with  loss  of  significance 
between  nerrly  equal  quantities  and  division  by  taro  (within  the 
computer's  finite  resolution;. 

Equation  (5.6)  it  derived  using  the  thermodynamic  relationship 


l£-  8  i£  -  Pl 

3v  38  P' 


in  conjunction  with  equation  (5.1).  A  more  direct  form  would  be 

(»,„) .  -x*x.  ,‘8j  htkllx  + 

Av  v",p/  lnvj+1  -  lnvj  2  *♦’ 


where 


A(lnv)  1 
— *7 — u  <v  “  ■  p, 
Av  V 


(5.17) 


(5.18) 


EJ  8J,k 


+  (8  -  8  v)  fl'  8j,k 
J,k  J,k+1  J,k 


(5.19) 


(5.70) 


S 


3 


Cfi  “  Cfi  +  <®  "  81+1  £>  fl 

8j+l  8J+l,i  J  1,1  8 


°S+ltt+l  »j+l. 
J+l,£+l  "  #j+M 


(3.21) 


The  cwo  form*  give  different  results  and  the  former  waa  choaan  only 
bacauaa  It  proved  to  be  cheaper. 

The  total  preaaure  la  the  sum  of  the  gwa  preaaure  given  by  aquation 
(5.1;,  the  radiation  preaaure, 

Pr  "  3  (3.22) 

and  the  pseudo- viscoue  pressure,  q,  which  la  an  artificial  aid  used  to 
fit  the  shock  In  a  smooth  manner  over  several  sonea  [3],  Two  forma  of 
pseudo-viscous  pressure  are  provided,  the  choice  of  vhfch  la  an  option 
provided  to  the  user.  The  form  linear  in  the  velocity  gradient  la 


k  n+1 
qi+Jj 


’l  pi+*j 


•  n+ij  -n+»j 

f 

*1+1  *1 

Ln+l  I 

\xi+l  / 

for  <  G, 


(5.23) 


for  AvJ^J  *  0, 


and  the  quadratic  form  la 


k  n+1 
qi+»i 


n  o°+1 

n2  Pi+Jj 


/  xn+1 \“  1  2 
;n-Mj  _  .n+Jj/  _ \ 

1+1  1  n+1 

L  '  *1+1  /  J 


for  AvlJ^  *  °* 


for  Av^  >  0, 


(5.24) 


Both  and  n2  be  altered  by  the  user  at  any  time  through  the  keyboard 

and  are  Initially  set  at  0.8  and  1.0  respectively.  The  parameter  o*0,l,2 
for  plane,  cylindrical  and  spherical  geometry  respectively.  The  sonic 
Vwljcity  la  given  b*-' 


k-1  n+1  .  A/,  fk-1  n+1  n+1  "1  ** 

4/3 !  Pi+*  Vl+^  I  * 


(5.25) 
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Note  chat  the  sonic  velocity  ia  a  function  of  the  total  pressure  and 
is,  therefore,  one  iteration  behind.  In  summary,  the  total  pressure 
is  given  by  the  sum 
k/  I k  /  .  n+Jj 


P(Q.p) 

li+h 


+  tj 


(5.26) 


5.2  Mean  Freo  Path 

Corresponding  to  a  conductivity  coefficient  in  thermal  diffusion 
programs,  a  dimensionally  equivalent  coefficient,  y  T,  is  calculated 
for  use  in  the  radiation  diffusion  term  of  equation  (2.6),  The  mean 
free  path,  T,  is  a  function  of  both  temperature  and  density  as  were 
the  preceding  energy  and  pressure  terms.  It  is  calculated  using  the 


relation 

r’£ 


(5.27) 


where  <,  the  opacity,  is  typically  calculated  through  the  interpolation 
of  tables  given  in  terms  of  density  and  temperature.  The  schema  is 
further  complicated  due  to  the  requirement  that  we  must  obtain  a  value 
for  A  at  the  interface  while  the  temperature  and  density  values  ere 
svallable  only  at  the  zone  mid-points.  Many  schemes  have  been  tried 
end  tested  [5,  6]  in  a  number  of  similar'  computer  programs. 

The  scheme  implemented  here  expands  upon  ideas  developed  in  the 
FF  program  15].  Thj  central  idea  is  to  try  to  get  e  better  estimate 
of  local  conditions  at  the  interface  and  then  calculate  a  mean  free 
path  based  upon  these  conditions. 

Of  prime  importance  was  an  interface  temperature  since  the  density 
changes  from  one  time  step  to  the  next  tend  to  be  overshadowed  by  much 


greater  changes  in  temperature.  Also,  the  opacity  tends  to  be  e  much 
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less  sensitive  function  of  density.  Thus,  density  .effects  were  some¬ 
what  Ignored  with  emphasis  upon  arriving  at  a  suitable  Interface 
temperature . 

Note  that  the  radiation  energy  aquation 
Er  -  ae4V  (5.28) 

A 

implies  that  <f>  -  0  is  a  measure  of  the  energy  In  a  zone  (as  far  as 
radiation  transport  is  concerned)-.'-  This  assumption-  is  valid  only  if  <p 
represents  the  whole  zone..  The  temperature  point  is; therefore,  located 
at  the  mass  center  rather  than  the- spatial' center  Of  the  zone.  This 
also  agrees  with  the  hydrodynamic-  differencing  as "explained  in  section 
2.2  and  given  by  equation  (2.7). 

The  radiation  energy  in  the -volume  between' mass' centers  about  a 
typical  Interface  is  simply 


a 

2 


+  (*V) 


(5.29) 


Also, 


2r  “  a(*V) 
i 

thus, 


♦i  " 


V.’a  +  v- 


(5.30) 


(5.31) 


This  temperature  is  exact.  Its  position  is. .unknown  though  and  one 
can  only  say  that  it  is  the  best-  estimate  which  can  be  obtained  in  the 
neighborhood  of  Interface.  In  the  end,  its' position  is  unimportant. 

What  is  important  is  that  it  provides  a  means'  of  obtaining  the  correct 
flux  across  the  interface  agreeing  with  empirical  cud  analytical  results. 
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Given  the  interface  temperature,  an  bpatity  is  calculated  for  each 


zone  in  the  neighborhood  of  the  interface  which  is  reflective  of  the 

respective  material  type,  temperature  and  density.  These  values  are 
+  — 

referred  to  as  ie^  and  where 

Ki  “  Kl9i»  Pi-fV  (5.32) 

and 


Ki  "  K(ei*  pl-lS)  (5.33) 

The  calculational  details  of  the  opacity  functions  are  given  in  the 
next  section. 


Note  that  the  opacity  is  a  measure  of  the  average  cross  sectional 


area  as  seen  by  a  photon.  Its  units  of  measure  are  area  per  gram  of 
material.  The  total  opacity  is,  therefore,  «cm,  where  m  is  the  ass 
of  the  material  which  has  opacity  tc .  Thus  in  our  case 


(Km)±  =« 


Vi-i,  +  K  t"m, 


and  since 


then 


(5.34) 


(5.35) 


-  V*  +  Vi** 

Kimi-Ss  +  Kimi-Hs 


(5.36) 


Since  each  zone  is  assumed  to  be  homogeneous  in  density,  the  mass 

center  is  also  the  volume  center.  Thus,  the  volume  between  x,  .  and 

i“Q 

xi  is  simply  .  This  fact  is  also  used  in  equation  (2,7). 
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An  average  flux  is  then  given  by 


F.  A 

i  3  i 


7  !li^) 


i4?2 


Xi 


(5.37) 


The  spatial  position  of  this  flux  level  is  not  fixed  and  is  not  known. 

It  becomes  important  only  when  the  radii  varies  between  x  ,  and  x 

i~=S  i+Jj* 

This  flux  value  is  also  valid  for  but  a  short  time.  Thus  care  and 
consideration  must  be  given  to  the  proper  selection  of  A,  the  area, 
and  at,  the  time  step  size,  which  are  used  to  give 


AQi  =  FiAiAt* 


(5.38) 


the  net  amount  of  energy  transported  from  one  zone  to  the  next  by 
radiation.  In  this  case,  the  area  used  is  the  actual  area  at  x^ 

It  will  in  general,  vary  as  does  y.±,  and  the  details  of  its  calcula¬ 
tion  and  those  of  the  zone  volumes  are  given  in  chapter  3.  The  time 
step  size,  At,  depends  upon  several  constraints  as  have  been  discussed 
xn  chapter  4.  Not  only  is  it  used  to  control  the  truncation  error  in 
—  and  consequently  —  ,  but  it  is  also  used  to.  limit  the  time  for 
which  a  calculated  flux  value  must  be  used.  Note  from  equation  (2.21), 
that  depending  upon  the  value  for  a,  the  new  and  old  flux  levels  are 
averaged  over  the  time  step. 

The  mean  free  path  at  the  boundaries  is  calculated  in  much  the 
same  manner  as  any  other  zone.  First  a  boundary  interface  temperature 
is  found  upon  which  an  opacity  is  calculated.  In  ..reference  to  equation 
(2.31),  this  boundary  temperature  is 


<i>0  =  <fi 


\  if  B0  *  0 


(5.39) 
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(5.47) 


or 
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\ (1  +  W 


a  +  yy2  + 

Here,  the  geometric  mean  free  path 

Ag  -  d. 


(5.48) 


(5.49) 


the  diameter  of  the  pipe  at  the  interface  in  question  and  Am  is  the 
so  called  material  mean  free  path  as  calculated  above  in  equations 
(5.36),  (5.42)  and  (5.46). 


5.3  Opacity 

The  opacity  tables  are  of  the  same  form  as  the  equation  of  state 
tables.  For  each  tabular  density  in  the  form  lWj ,  there  is  a  sequence 
of  temperature  and  opacity  pairs,  (lnS,  ln<)  such  that 

J 


lnvx  <  lnv2  <  ...  <  lnv  , 


ln9j>i  <  ineJj2  <  ...  <  ln9j>K  (5.S1) 

with  J,  K  >  1.  These  tabl.s  are  also  generated  by  other  programs  [8], 

and  a  graat  number  are  available  corresponding  to  the  aquation  of  stat. 
tables. 


The  opacity  table  is  interpolated  for  Ins  from  which  the  mean  fr.e 

path  is  calculated  in  equation  (5.36).  The  interpolation  formula  is 
simply 

lnv  -  lnv 

ln«(6,p)  inv.j+1  -  lnvj  ^lnKj+l  “  lnKj)  +  , 


(5.52) 


so 


where 


ln0  -  ln0 


ln< 


iiiL 


J+1  ln6j,k+l  '  ln6j,k  (  lnKJ*k+1  ~  lntCJ.k)  +  lnKJ  ,k 


(5.53) 


ln0  -  ln0 


lnvi '  <invi.m  -  +  iru'j+i,i<5'54) 


assuming  that 

lnVj  £  lnv  s  with  1  s  j  <  J, 

ln0j  k  S  Jn0  S  ln0j  k+1  with  1  *  k  <  K* 
ln0j+l,i  *  ln0  S  ln0j+l,£+l  with  1  *  A  <  K‘ 


(5.55) 

(5.56) 

(5.57) 


As  was  the  case  with  the  equation  sf  state  tables,  the  equations  are 
also  used  as  extrapolation  formulae  if  the  density  and  temperature 
values  are  outside  of  the  table  limits.  This  does  cause  some  undeslr** 
able  problems  especially  in  the  low  temperature  regions. 
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6. 

SOURCE  TERMS 


In  addition  to  the  boundary  pressure  profiles  *hich  havs  already 
been  discussed,  the  user  may  impose  time  dependent  temperature  or 
energy  profiles  upi-n  any  set  of  ;.ones  within  the  problem.  Both  may 
not  be  done  simultaneously,  but  the  user  may  change  from  one  to  the 
other  dynamically  as  the  solution  progresses  dependent  only  upon  his 
problem  restraints  and  needs.  He  may  also  modify  or  otherwise  change 
the  one  he  is  currently  using.  Both  require  time  and  space  dependent 
profiles.  The  space  dependence  is  on  a  zone  basis  and  it  is  very 
easily  changed.  It  has  been  referred  to  previously  as  the  set  of 

power  factors  P,  ,  i  •  0,  . . . ,  N  -  1. 

14*5 

The  time  dependent  temperature  or  energy  profiles  are  assumed  to 
be  in  tabular  form,  but  functional  forms  are  simply  provided  for  by 
the  user  supplying  his  own  subroutines. 


6.1  Energy-Time  Profile 


An  energy  table  of  energy-time  pairs  in  the  form  (E^,tj)  may  be 
specified  by  the  user  when  the  program  is  started,  or  he  may  enter  or 
alter  the  table  dynamically  during  program  execution.  He  need  only 
be  careful  that 

(6.1) 

where  t"  is  the  current  problem  time,  and  that 
•  •  •  t  . 


s  tn, 
.  n 


C1S  C2 


(6.2) 
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with  J  &  20.  The  source  energy  for  each  zone  la  then  given  by 


ASn+!*  -  P 
i-Hs  f 


(E 


-  EJ 


i±i _ !il  Atn+*s 


i-Hs  (tJ+l  "  tj) 


(6.3) 


where 


v  '■ « ‘n+l  ‘  Vi  • 


(6.4) 


If  the  time  step  size  is  so  lcrge  as  to  extend  over  parte  of 
more  than  one  table  interval,  then  the  source  integration  ia  done 
piece-wise  in  order  that  the  problem  will  reflect  the  correct  total 


energy  at  all  times.  If 
<n+1 


(6.5) 


then  the  table  is  simplt  extrapolated  in  a  linear  fashion  using  the 
last  two  entries.  This  of  course  implies  that  the  table  must  have  a 
minimum  of  two  entries.  Note  that  the  table  need  not  even  be  piece- 
wise  continuous  and  that  multiple  entries  may  be  given  for  the  saaie 
time.  This  facilitates  step  as  well  aa  ramp  energy  excursions.  For 
example,  if 

n+l 


Vi 1  ‘" < 


then 


Vi <  c 


£  t 


j+2  • 


(6.6) 


C  ■  '« 


<«j  -  f>  ♦  <Vl  - .  > 


LCJ  ~  Cj-1 


+  (tn+i  - 1, . , ) 


Cj+2  "  CJ+.l 


J+1‘ 


(6,7) 


Care  should  be  taken  to  insure  tl  t.  j*  t_. 

J  J 
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6.2  Temper \ tura-Time  Profile 

The  stricture  of  a  temperature  profile  is  not  unlike  that  which 
has  been  giv  in  for  an  energy  profile.  It  is  a  profile  of  temperature 
time  pairs  t,l  the  form  (0^,  tj)  specified  in  the  same  manner  as  is 
the  energy  profile.  The  purpose  of  this  profile  is  to  impose  a  temp¬ 
erature  over  a  set  of  zones  within  the  problem.  These  zones  thus 
act  much  like  a  time  dependent  heat  reservoir,  be  it  a  sink  or  a 
source.  The  imposition  of  the  profile  over  a  specific-  set  of  zones 
is  done  through  the  power  factors.  If  P  -  1,  then  this  is  used 

as  a  flag  to  indicate  that  the  temperature  of  this  zone  is  determined 
by  the  profile.  Likewise,  if  P  ■  0,  then  the  temperature  of  this 

zone  is  not  dependent  directly  upon  the  profile. 

Rather  than  simply  "overloading"  the  specified  zones  with  the 
temperature  as  given  by  the  profile,  it  is  desirable  to  calculate  the 
source  term  necessary  to  give  this  temperature  for  each  zone  in  ques¬ 
tion.  It  is  then  possible  for  the  user  to  monitor  the  amount  of 
energy  being  "dumped"  into  (or  out  of)  the  problem.  He  is  also  able 
to  monitor  the  overall  problem  energy  balance. 

The  first  step  toward  calculating  this  source  term  is  to  inter¬ 
polate  the  profile  for  the  current  temperature: 


^J±L^rL  <‘“+1  -  v  * ... 

j+i  i  j 


where 


n+1 


£  t 


(6.8) 


j+1  * 


(6.9) 
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Then,  the  extrapolation  procedure  depicted  in  figure  3  and  given  by 
equation  (2.53)  ie  modified  as  follows: 


0  n+Jl 

°i+Jj 


v^(6?„  -6^>+9^  ifPf 


At 
0  if  P 


n-S  vwi^ 

i  0. 


i-ffl 


“  0, 


■i-Hs 


(6.10) 


After  obtaining  an  estimate  of  the  zone  temperatures  as  indicated 
above,  the  temperature  dependent  terms  are  updated  as  indicated  by 

equations  (2.54),  (2.55),  (2.56),  (2.57),  (2,58)  and  (2.59).  The 

k 

source  tern  ASi+^  is  then  calculated  (and  recalculated  for  each  iter¬ 
ation)  through  equations  (2.36),  (2.37)  and  (2.38)  by  replacing  each 


kAn+l  . 
*i+4  b* 


e4  if  pf  +  o. 
*1+4 


Note  also,  that  source  terms  are  cal¬ 


culated  only  for  those  zones  for  which  P,  i  0.  In  actual  practice 

‘i+Jj  * 

Pf  is  used  as  a  multiplicative  factor  whan  solving  for  from 

i+*l  i+*i 

each  equation. 

c  >  then  equation  (6.8)  is  used  ns  an  extrapolation 


formula. 


ss 


7. 

COMPUTER  CALCULATIONS 

A  number  of  calculations  wore  done  to  illustrate  the  effect  and 
tlie  worth  of  various  program  parameters  and  to  also  validate  the  method 
of  soluti°n.  In  addition,  F3,  a  one-dinansional  program  at  Los  Alamos 
[5],  was  used  to  do  some  of  the  same  cnlculatir  a  three  way  compari¬ 
son  was  then  made  with  the  analytical  solutions. 

The  program  developed  from  the  method  given  i**  mis  paper  is 
referred  to  as  HYRAD1.  It  has  been  programmed  to  run  on  the  PDP-10 
-omputer  at  the  University  of  Utah  via  the  remote  terminal  at  Montana 
.State  University  and  also  on  bath  the  CDC-6600  and  CDC*.’<iOO  computers 
at  Los  Alamos.  The  6600  und  7600  versions  are  identical,  being  FORTRAN 

program-,,  but  are  restricted  at  this  time  to  plane  geometry  problems 
only.  .  . 

7.1  Radiation  Diffusion  Calculations 

A  series  of  calculations  is  presented  here  which  was  used  to 
validate  the  radiation  diffusion  calculations  and  to  illustrate  the 
effect  and  worth  of  several  program  parameters.  A  complete  description 
of  the  problem,  often  referred  to  as  the  Marshak  Wave  Problem,  may  be 
found  in  .9].  The  basic  elements  are  as  follows: 


1. 

Constant  specific  heat. 

l£ 

aa  * 

2. 

No  energy  in  the  radiation  field. 

3. 

Constant  density,  p. 

4. 

An  opacity  *c  ■  Kooa0~fi 

5. 

Constant  driving  temperature,  0^, 

6. 

Plane  geometry. 
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Thus,  for  simplicity: 


Of 

00 


nc 

3 


_4_ 

n+4 


<■  ■  ’O 

'Q  U 

a  -  l 
60 

*?+«,  ■ 10'6 ' 0 


1 


(7.1) 


In  addition,  a  unit  cross  sectional  area  was  assumed.  Several 
values  for  8  were  tried  in  addition  to  varying  length  zones,  etc. 
Corresponding  to  the  tables  published  in  [9],  the  space  and  tempera¬ 
ture  values  were  normalized  as  follows: 


(7.2) 


c“,  .34 

**  V7 

'rt  •  “rt 1  %  -  •  <7-3> 

The  ralculaclons  could  then  be  checked  at  any  tine  for  which  there  was 
n  xi+*j  whlch  8av«  rise  to  a  6^  that  appeared  in  the  published  tables. 
This  was  not  difficult  since  tabular  values  were  given  in  increments 
of  0.05  for  £  starting  with  C  ■  0  and  t  -  0q  /  0^  ■  1  out  to  where 
t  -  0.  No  interpolation  was  necessary  nor  desired. 

The  purpose  of  these  calculations  is  to  follow  the  diffusion  of 
n  radiation  wave  driven  by  a  constant  boundary  temperature.  Its 
progress  is  then  checked  at  various  times  with  respect  to  its  position 
and  shape.  Figures  4  through  10  give  the  position  of  the  wave  at 
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di/.ferent  times  as  a  temperature  prof  11a.  Zn  theta  display  conaola 
pictures,  the  temperature,  prassura  and  danslty  lncraasa  to  rha 
right,  and  tha  x  coordlnata  lncraasas  In  the  vartlcal  dlractlon. 

Tha  density  curve  at  tha  right  la  constant  over  all  time  and  tha 
horizontal  lines  are  Indicative  of  the  zone  boundaries. 


Figure  5 
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Figure  7 
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Note  that  a  constant  density  prohibits  any  hydrodynamic  activity 

and  that  this  problem  is  purely  a  radiation  diffusion  calculation.  The 

first  serids  of  problems  was  done  with  6-0  and  the  zone  width  Ax^ 

.2  cm.  Thus,  X  ■  “  -  .1  cm,  only  half  the  zone  width.  The  initial 

0 

temperature  distribution,  0°^-  10~6,  is  sufficiently  close  to  zero 

and  is  very  near  the  calculational  limit  of  the  PDP-10  for  $  -  e4  =  io-24 
and  |~  ^  0.2  —  0.274442. 

Table  1  below  gives  the  exact  solution  and  corresponding  values 
for  t  as  calculated  by  HYRAPl  and  F  at  common  values  of  £  for  t  =  36  sh. 
The  tabular  entries  are  rounded  to  five  digits,  and  it  is  important  to 
point  out  that  even  in  double  precision,  some  significance  in  the  fifth 
digit  is  the  most  one  can  expect  on  the  PDP-10  computer.  For  HYRAD1, 
a  -  ]/2,  v  -  1,  e  -  10-5,  Tfac  -  .10,  and  At0  =  10~7  sh.  While  for  F3, 
Tfac  =  .10  and  At0  =  IO-8  sh. 
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c 

T 

HYRADl 

F3 

.05 

.98803 

.98803 

.98803 

.15 

.96273 

.96272 

.96272 

.25 

.93538 

.93537 

.93537 

.35 

.90563 

.90562 

.90562 

.45 

.87304 

.87303 

.87303 

.55 

.83699 

.83697 

.83697 

.65 

.79661 

.79658 

.79658 

.75 

. 75061 

.75056 

.75056 

.85 

.69692 

.69685 

.69685 

.95 

.63187 

.63177 

.63177 

1.05 

.54763 

.54754 

.54755 

1.15 

.42047 

.42008 

.42009 

1.25 

0  .  , 

.00638 

.00637 

S.E.(j) 

1.5490 

1.5489 

1.5504 

I.E.(j) 

1.5490 

1.5489 

1.5489 

#  time  steps 

2316 

*327 

avg.  #  iterations  , 

6.1 

1 

calculation  tipe  (min) 

* 

@  49  sh 

7600 

1 

.4 

PDP-10 

* 

64 

Table  1:  >  HYRADl  and  F3  Calculations  As  Compared 
With  The  Analytical  Solution  For  K  ■  10. 

-  I  '  ‘  *  ‘  ' 

The  next  thrfee  Sets  of  tables  and  figures  illustrate  the  difference 

<  \ 

caused  by  variations  in  o,  v,  Tfac  and  €  as  calculated  by  HYRAbl.  In 
table  2  and  figure  11  the  only  changes  in  t  occur  between  o  ■  h  and 
a  ■  1  for  6  »  10'  and  T-  .  ■  .10.  For  a  ■  0,  the  calculations  go 

XAC 

unstable  after  6  shakes  when  At  >  1.5  »  10  The  results  at  4  shakes 

are  impressive  though,  being  much  better  than  those  for  a  "*3»  The 

severe  time  step  tilze  restriction  for  this  parameter  value  is  less 
than  desirable  and  no  further  calculations  were  tried  at  this  time 
with  o  -  0.  The  effect  of  the  extrapolation  parameter  v  is  to  reduce 
the  number  of  iterations  required  for  each  time  step  advancement'.  For 

"■  *  ‘  ,  ’  *  *  .  ,  i  * 

this  and  most  problems  it  is  most  effective  at  unity.  For  some  hr°hlems, 

adjusting  v  can  reduce  the  number  of  iterations  by  as  much  ah  a  third. 
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Running  tlmaa  at  glvtn  for  tht  POP-10  art  approxlmata  and  vary  aavaral 

minutes  depending  upon  tha  machine  activity  (i.e.  tht  number  of  tlmaa 

tht  program  mutt  bt  swapped). 

Number  of  Avg.  Numbtr  of  PDP-10  Computtr 

a  v  Tima  Ste.pt  Itaratlona  •'  Tima  (Min) 

0  1  (unstabla) 

>5  1  2316  6.1  64 

1  1  2304  6.1  64 

1  .5  2304  7.1  74 

1  0  2304  7.3  75 

Tabla  2:  Variations  in  tha  Tima  Diff trancing  and 
Extrapolating  Paramaters,  o  and  v,  at  t  -  49  ah. 
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Table  3  and  figure  12  show  the  relative  insensitivity  of  result* 
upon  the  tine  atep  *lze  restriction  factor  x£||-  for  o  -  1,  u  -  1  end 
€  "  10  5*  Thi*  al#0  holds  true  for  a  -  %.■  Note,  however,  that  it 
daes  dramatically  effect  the  calculation  tine  directly  by  Halting  the 
tine  step  slse. 


T, 

fee 

Number  of 
Time  Steps 

Avg.  Number  of 
Iterations 

Avg.  At 
(•h) 

PDP-10  Compute?: 
Time  (Mii) 

.01 

27993 

2.2 

.00175 

295 

.05 

5478 

4.2 

.00895 

100 

.10 

2661 

6.1 

.01841 

70 

.15 

1720 

7.8 

.02849 

52 

.20 

1248 

9.6 

.03926 

46 

.25 

966 

11.4 

.05072 

42 

.  30 

774 

13.2 

.06331 

41 

.35 

641 

15.1 

.07644 

38 

.40 

540 

16.9 

.09074 

36 

Table  3i  Variation*  in  the  Temperature  Restriction 
Factor i  Tfcfef  at  t  -  49  sh. 


T  HI  .Rip 
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Table  4  and  figur.  „  lllu.tr.,..  th.  .ignite.  #f  t„. 

..nc  pn«.t«  .nd  tfM.  oa.  vauM  of  €  i<m 

'  an  or  agul  t0  10  ,  th.  raaul„  .r.  un.ff.cted,  Thi.  i. 

.l»o.  4-5  df,lt.  of  elgniflcenca  i.  .u  «b.t  c.„  b.  obt.i„.d  on  th. 
PDP-10.  What  is  lnpr.s.lvo  1,  th.  reasonably  good  r.«ult.  for  €  -  Hf1. 

Tak*  not*  also  of  the  variant. 

“C*  in  the  average  nunbar  of  iteration,  par 

atop  and  it.  affect  upon  the  calculation.!  tina. 


10"7 

10-6 

10-5 

10-4 

lO"3 

10-2 

io-i 


I.E. 

(J) 


1.806998 

1.806998 

1.806998 

1.806997 

1.806966 

1.806708 

1.806070 


S.E. 

(J) 


1.806998 

1.806998 

1.806999 

1.806999 

I. 207028 

J . 807243 
1.807762 


Numbar  of  Average 
Time  Steps  Numbar  of 


Iterations 


2662 

2662 

2661 

2661 

2661 

2662 

2671 


8.8 

7.4 

6.1 

4.8 
3.2 

1.8 
l.C 


PDP-10 
Computer 
Tiaa  (Min) 

90 

80 

70 

58 

39 

28 

23 


actual  1.807152  1.807152 

Tabla  4;  Variations  in  the  Convergence  Parameter,  €  ,  at  t  -  49  th. 
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Several  other  calculations  war*  dona  with  other  values  of  8. 
labia  5  gives  the  results  for  a  calculation  at  36  ah.  don*  with  8-3. 

3 

Far  this  case.  Cy  -  ~  ^  -  .156825,  X  -  a  -  h,  v  -  1,  €  -  10*5, 
Xfau  "  *10  *nd  At°  "  10  7  8h*»  and  for  p3»  Tfac  “  -iO  “d  6t°  -  10'8  ah. 


T 

HYRAD1 

F3 

.05 

.94296 

.99297 

.99297 

.15 

.97792 

.97793 

.97793 

.25 

.96139 

.96141 

.96141 

.35 

.94304 

.94307 

.94307 

.45 

.92242 

.92246 

.92246 

.55 

.89887 

.89891 

.89892 

.65 

.87136 

.37141 

.87142 

.75 

.83820 

.83828 

.83828 

.83 

.79616 

.79635 

.79634 

.91 

.73785 

.73849 

.73845 

1.05 

.63704 

.63945 

.63933 

1.15 

0 

.00021 

.00015 

S.E.(J) 

.$9914 

.90917 

.90989 

I.E.(j) 

.90914 

.90917 

.90887 

#  time 

steps 

2457 

2468 

avg.  # 

iterations 

5.9 

1 

calculation  time 

@  49  sh  CDC-7600 

L  PDP-10 

1  min 

68  min 

•  4  ciu 

Table  5:  HYRAD1  and  F3  Calculation*  as  Compared  with 
the  Analytical  Solution  for  k  -  10/0^. 


Another  series  of  problems  were  run  for  a  variety  of  zone  widths 
from  0.1  cm  to  200  cm.  The  difference  in  the  results  wer*  significant 
and  Improved  with  smaller  Ax's  approaching  the  mean  free  path.  This 
just  confirms  one's  Intuition  that  the  linear  space  derivative  approxi¬ 
mation  for  does  not  adequately  describe  the  flux  terms  for  coarse 
zoning  even  for  materials  which  exhibit  constant  or  nearly  constant 
opacities.  So  much  concern  and  work  has  been  directed  toward  calculating 
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*  mean  fraa  path  at  tha  intarfaca  batween  zones,  whereas  it  should 
probably  ba  directed  at  arriving  at  a  batter  flux  term  as  a  whole. 

Of  course,  thia  is  difficult  to  do  properly,  for  there  is  just  not 
enough  information.  In  most  cases,  it'  is  usually  better  to  zone 
finer.  No  calculations  ware  done  with  zone  widths  less  than  a  mean 
fraa  path. 

One  interesting  result  obtained  from  the  variable  zone  width 
calculations  was  tha  invariance  of  the  results  on  a  time  step  by 
time  step  basis.  More  precisely,  if  the  calculation  time  is  nor¬ 
malized  with  respect  to  the  square  of  the  multiplicative  difference 
in  zona  widths,  then  the  results  at  comparable  normalized  times  are 
the  same.  For  example,  the  temperature  profile  at  t  -  9  sh  for 
Ax  -  0.1  ia  the  same  as  for  t  «  36  sh  for  Ax  -  0.2,  but  the  latter 
calculation  includes  exactly  twice  as  much  energy. 

7.2  Hydrodynamic  Calculations 

In  this  section,  a  series  of<  hydrodynamic  calculations  are  pre¬ 
sented.  Two  different  problems  were  Investigated.  .  The  first  is  h 
shock-tuba  problem  which  is  described  in  [10]  and  [11] .  This  calcula¬ 
tion  follows  the  shock  formed  by  a  high  pressure  gas  expanding  into  a 
low  pressure  gas  aonfinad  in  a  long  small  radius  pipe  or  tube.  The 
problem  assumes  an  ideal  gas  for  which 
e  -  p/p(A-l) 

P  -  PR6, 


with 


(7.4) 

(7.5) 


R  -  C  -  C  . 
P  v 


P  "  v*  (7.7) 

9  9 

For  air,  A  -  1,4  and  R  -  286.793  m  /sec2  °K  assuming  29  grams/mole 


of  air. 


From  equations  (7.4)  and  (7.5), 


R/A-l, 


716.983  m2/sec2  °K, 


(7.8) 


3v  “  u*  (7.9) 

i  i  e 

For  the  one  region  of  gas,  p  -  666.447  kg/m-sec2  and  p  -  .0077459716 
3 

kg/m  ,  and  for  the  other  region  p  -  1.1823  x  107  kg/m-sec2  and 

3  ■'  ^  •-  ii 

p  ■  137.413  kg/m  .  The  initial  tfemperature  for  both  regions  is  300°K. 

The  first  Calculation  done  was  with  a  constant  Ax  ■  .0254  m, 

-6 

At  -  1.25  x  10  see,  r  «*  .0254  m  and  an  artificial  viscosity,  factor 
“  0.8.  The  results  of  thijs  calculation  match  the  results  for  the 
K0  and  PUFL  programs  as  reported  in  [10] .  The  next  calculation  was 
done  with  n2  “  2.  Again  the  overall  result  was  the  same.  However, 
when  the  region  at  and  behind  the  shock  front  Was  examined  in  detail , 

t  •  »• 

the  quadratic  form  of  artificial  viscous  pressure  exhibited  rather 
large  oscillations  as  displayed  in  the  pressure  profile  given  in 
figure  14.  Notice  in  contrast  the  smooth  pressure  profile  for  the 
linear  form  of  the  viscous  pressure  term,  A  combination  Was  then 
tried  with  *  1  and  n2  -  3  as  suggested  in  [lb].  These  results 
were  identical  with  the  linear  case  except  right  at  the  shock  front  l"| 
which  is  not  as  sharp  and  is  extended  over  several  more  zones. 


Another  series  of  calculations  were  done  with  a  variable  At.  j?or 
this  series,  the  major  time  step  restriction  from  thope  discussed  in 
chapter  4  is  the  hydrodynamic  restriction  determined  by  the  parameter 
Hf ac ’  Two  calculations  were  done  with  Hfac  -  .05  and  .10  and  there  was 
little  significant  difference  in  the  results  as  shown  in  figure  15. 
However,  there  was  a  great  difference  in  the  calculational  efficiency 
as  indicated  by  table  6. 
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Constant  At 


avg.  time  step  size  (sec)  1.25  x  iq"*6 

total  number  of  time  steps  3120 

avg.  number  rf  iterations 
per  time  step  4.7 

PDP-10  calculations!  time  (min)  228 


Hfac  “  *05  H^ac  ■  .10 

1.915  x  10"6  4.004  x  io“6 
2037  974 


5.7  7.5 

151  92 


K.E.  (j) 
I.E.  (j) 


2.02030x10^ 

-2.02063x10^ 


2.01989x10^ 

-2.02033X104 


2.02013x10^ 

-2,02037x10^ 


Table  6:  Relative  Calculational  Efficiency  of  Constant  vs.  Dynamic  Time 

Step  Sizes  at  t  -  .0039  sec. 


Another  calculation  done  was  with  a  constant  At  but  with  Ax  ■  .0508 
meters,  double  that  of  the  previous  calculations.  As  was  the  case  with 
the  pure  radiation  diffusion  calculations,  the  solution  was  very  sensi¬ 
tive  to  the  zone  width.  Figure  16  gives  the  difference  in  the  pressure 
profiles  at  the  time,  t  -  .0039  sec. 
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—*wu*otion8 


radlus  tube.  Here 

wa®  Initially  .pacified  as  the  ’  C°”1C  8eCtl°n'  °r  ttUBtrm 

container  of  the  lnu  nv 

■-  M«h  pressure  f!as  was  then  u  "  8“' 

before.  In  thl  all°wed  t„  expand  into  this  w1m.  aa 

rate  one  Inch  ,  Se“1<m’  £h<  radlUS  lnCreMe<1  at  the 

The  calculations  started  with  init^i 

h  initial  one  inch  (.0254  ml 
«ne  widths  and  specific  a  ,  * 

P  cified  a  constant  At  -  1.25  x  i0~6  8ec. 

The  first  problem  was  run  -u 

-trrrirrir 

nr  zone  center.  ^  “lal  "“‘Point 


Flg-  171  Zone  "“‘Point  Definitions 


The  results  of  the  calculations  are  given  in  figure  18  at  t  ■  .0005 
and  .001  seconds  in  the  form  of  temperature  profiles.  The  temperature 
was  used  here  in  lieu  of  the  pressure  because  it  did  not  vary  over  as 
wide  a  range  of  values,  but  is  directly  proportional  and  indicative  of 
the  pressure.  In  both  cases,  figure  18  shows  that  the  peak  temperatures 
in  the  center  of  volume  problem  are  greater  at.d  are  advanced  further 
down  the  tube  than  those  for  the  axial  center  problem. 

Table  7  gives  the  program  statistics  for  this  set  of  calculations. 
Note  that  for  this  geometrical  configuration,  the  kinetic  energy  as 
given  by  equation  (2.64)  does  not  match  the  loss  in  internal  energy. 

This  was  not  the  case  for  the  constant  diameter  problems  as  shown  in 
table  6.  At  about  .00075  sec,  the  problem  became  Courant  limited 
(see  section  4.1)  because  the  shock  was  encountering  larger  and  more 
massive  zones  as  the  pipe  diameter  increased.  Thus,  the  zones  tended 
to  pile  up  into  smaller  widths. 


CiS 


Axial  Center 

Volume  Center 

time  step  size  (sec.)  , 

@  .0005  sec. 

1.25  x 

1°J 

1.25  x  io“; 

@  .001  sec. 

6.48  x 

10  7 

6.11  x  10"' 

total  number  of  time  steps 

@  .0005  sec. 

401 

401 

@  .001  BtiC. 

9D8 

933 

avg.  number  of  iterations 

@  .0005  sec. 

2.7 

2.7 

@  .001  nec. 

2.9 

2.9 

PDP-10  calculational  time  (min) 

@  .0005  sec. 

9.9 

10 

@  .  001  sec . 

23.3 

23.6 

K.E. 

@  .0005  sec. 

3.23157 

x  10, 

3.19644  x  10, 
7.11488  x  lbJ 

@  .001  sec. 

7.22040 

x  10J 

I.E. 

@  .0005  sec. 

-3.14898 

x  10* 
x  10* 

-3.11402  x  10 * 
-6.76969  x  loJ 

@  . 001  sec . 

-6.82762 

Table  7:  Program  Statistics  for  Axial  vs.  Volume  Zone  Centers 


The  second  problem  investigated  in  the  hydrodynamic  series  is 
known  as  the  Von  Neumann  point  source  problem  [12],  This  purpose  of 

this  problem,  also  known  as  the  blast  problem,  is  to  calculate  as  a 

» 

function  of  time  the  blast  radius  propagating  from  the  blast  point  in 
spherical  geometry.  The  Von  Neumann  solution  to  the  blast  problem 
is  given  in  [12]  and  numerical  solutions  are  given  in, [13]  and  [14] . 

This  problem  is  similar  to  the  shock  tube  problem  in  that  it  assumes 

:  .  :  * 

an  ideal  gas  equation  of  state  (equations  7.4  through  7.9),  but  differs 
in  that  a  large  amount  of  energy  is  released  at  (near)  the  center  of  a 
spherical  volume. 

Eilers  and  Whitfill  at  Los  Alamos  have  been  using  this  problem  to 
validate  numerical  integration  techniques  and  to  establish  parametric 
values.  In  particular,  they  provided  several  calculations  done  With 


the  F  program  for  comparison  with  HYRAD1  and  the  analytical  solution 
given  by  [12].  One  such  solution  is  given  in  figure  19  together  with 
that  from  HYRAD1  and  variations  therein.  Figure  19  is  a  plot  of  the 
absolute  difference  between  the  actual  blast  radius  and  the  calculated 
blast  radius  as  a  function  of  time.  For  this  series,  AxJ  -  30.5  cm, 

Y"1,2’  ¥e“5’  Tv  "  R  “  1  (meaning  6  absorbs  the  actual  value  of  R) 

eWs  ■  10"6'  "l  ■  -8>  rfac  ■  -10-  8p  '  •0;'1>  Hfac  -  .01,'  At0  -  10'7  and 

(the  energy-in  at  the  blast  point)  E.  -  41b 5  . 

in  j 

At  this  point,  it  is  important  to  point  out  the  difference  in 

the  acceleration  terms  as  calculated  in  F3  [5]  with  those  calculated 
3 

by  HYRAD1.  F  uses  the  following  difference  equation  to  arrive  at  a 
numerical  solution  to  the  momentum  equation: 


1  ^ml-k  +  m1+ 
which  is  equivalent  to 


i  +  «1+4s) 


I  "  *1+%  xi-h 

I4i.  “  **(*.  + 


R  calc  -  R  actual 


Fig.  19:  Comparison  of  hydrodynamic  differencing  schemes 
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Equation  (7.11)  is  now  in  the  same  form  as  given  for  HYRADl  in 
equation  (2.7).  The  quantity  x^  as  given  in  (7.12)  is  said  to  be 
located  at  the  radial  center  in  contrast  to  the  center  of  volume  as 
specified  for  HYRADl  by  equation  (3.28).  The  difficulty  with  equation 
(7.11)  is  the  location  of  the  volume  given  by  V^.  The  denominator 
term  Implies  that  it  should  include  half  the  mass  of  each  adjacent 
zone  if  a  uniform  density  is  assumed  as  is  done  in  HYRADl.  In  con¬ 
trast  then. 


V1  ■  4irxiAxi, 


(7.13) 


which  when  expanded  becomes: 


Vi  “  47rL*t(x1+is  +  2x1+iix1_^  +  (*14ij  - 

If  was  considered  to  be  the  volume  between  x.  and  x..  ,  then 


(7.14) 


^i  “  4irC1/3(*144s  +  x14isxi_la  +  “  Xi-Js^ 


(7.15) 


It  is  not  difficult  to  show  that 

^i  "  Vi  "  3  *Xi44s  “  xi-^3,  (7.16) 

A  calculation  done  with  V±  in  lieu  of  V±  is  given  in  figure  19  and  is 
labeled  AVG^, 

Another  differencing  scheme  was  suggested  by  L;  A.  Schmittroth, 
and  it  proceeds  as  follows  from  integrating  the  momentum  equation  by 
parts; 

■  _  i£ 

^Jt  3x 


(7.17) 


Sl  '  *'  +  *■£)  <Pl-li-W  <7i22> 

In  oontrast  to  equation  (7.10)  which  uses  the  area  at  the  interface  $ 
equation  (7.22)  uses  the  average  of  the  areas  at  the  mid-point  of  the 
aonee  adjacent  to  the  interface.  The  results  of  the  calculation  using 
difference  equation  (7.22)  is  given  in  figura  19  and  is  labeled  AVG^i. 
Similarly,  if  we  let 


*Xi+»S  "  Xi-i5)  ’ 


(7.23) 


then 


HYRAD1  was  then  modified  to  calculate  the  zone  centers  and  accel¬ 
eration  terms  as  given  by  equations  (7.10)  and  (7.12)  for  the  F3  program. 
The  results  are  also  shown  in  figure  19. 

7.3  Combined  Radiation  Diffusion  and  Hydrodynamics 

The  results  of  a  combined  radiation  diffusion  and  hydrodynamic 
problem  are  shown  here  in  figures  20  through  25.  This  problem  is  a 
multi-material  calculation  In  which  one  material  is  heated  and  then 
transfers  this  energy  to  the  other  by  radiation.  The  initial  trans¬ 
fer  of  energy  creates  a  shock  in  the  second  material  as  indicated  by 
the  formation  of  a  spike  in  the  pressure  profile  shown  in  figures 

20,  21,  and  22.  This  influx  of  energy  causes  the  Second  material  to 

'■  1 

expand  one  zone  at  a  time  back  Into  the  first  material  and  also  for¬ 
ward  with  the  shock  causing  an  increasing  number  of  zones  to  be  com¬ 
pressed  at  and  immediately  behind  the  leading  edgfe  of  the  shock  wave. 

Note  that  from  figure  23  cn,  the  shock  Wave  has  overtaken  and  is 
accelerating  ahead  of  the  diffusion  wave.  This  is  indicated  by  the 
formation  of  the  additional  step  Within  the  temperature  profile. 
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7.4  Summary 

The  selection  of  computer  calculations  included  in  this  section 
is  not  exhaustive  but  rather  illustrative.  A  great  deal  of  work 
remains  to  be  done,  particularly  within  the  areas  of  validating  the 
non-linear  approximation  techniques,  the  hydrodynamic  activity  for 
various  geometrical  shapes  and  multi-material  effects.  Detailed 
studies  also  need  to  be  made  into  the  effects  of  the  equation  of 
state  and  opacity  calculatlonal  schemes.  For  example,  figure  26 
illustrates  a  difficult  yet  convergent  time  step.  Notice  the  plot 
of  the  relative  convergence  error  at  each  iteration.  The  initial 
decay  is  a  rather  moderate  exponential  response  which  changes 
abruptly  after  the  eleventh  iteration. 

The  other  plots  on  this  figure  show  the  change  of  several  zona 
variables  over  the  course  of  the  time  interval.  The  zone  plotted  ie 
that  which  had  the  largest  relative  convergence  error  at  iteration  15, 
the  last  iteration.  The  values  plotted  (from  top  left)  are  ~, 

P*  A,  the  luminosity  (flux  times  the  area)  into  one  end  of  the 
soae  end  out  of  the  other,  and  a  somewhat  confusing  and  uninformative 
convergence  graph  which  ha**  since  been  dropped.  The  initial  vertical 
line  on  the  ~  plot  (labeled  DEDT)  le  simply  a  line  drawn  from  tho 
base  line  to  the  value  of  as  it  was  three  time  steps  ego.  The  next 
two  time  step  values  ere  then  plotted  with  e  line  between  them  giving 
e  short  curve  which  le  representative  of  the  behavior  of  the  derivative 
over  the  previous  three  time  steps.  Another  vertical  line  ia  than 
drawn  from  the  base  line  up  to  the  first  iterative  Value  of  the  term 
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for  the  current  time  step  and  subsequent  iterative  values  are  represent¬ 
ed  by  the  curve  following.  This  format  is  repeated  for  the  p  and 
X  terms  (labeled  DEDV,  P  and  LAMBDA  respectively).  The  curves  for  0 
and  the  luminosity  represent  only  the  iterative  values  for  the  current 
time  step  and  (unfortunately)  do  not  include  values  for  previous  time 
steps. 

Similarly,  figure  27  is  a  picture  of  a  time  step  which,  while 

I 

convergent,  displays  rather  erratic  behavior.  Note  that  it  initially 
starts  to  converge,  abruptly  diverges  and  then  converges.  Notice 
tha  corresponding  graphs  of  the  behavior  of  various  zone  quantifies. 
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Figure  28  gives  a  snapshot  of  a  time  step  which  doesn't  converge 
after  thirty  iterations.  It  also  looks  as  if  it  will  never  converge 
and  appears  to  oscillate  every  five  iterations.  Notice  the  insensi¬ 
tivity  of  the  mean  free  path  and  one  of  the  luminosity  terms.  The 
other  luminosity  term  is  somewhit  affected,  and  the  pressure  even  less 
affected  by  the  fluctuations.  Clearly,  it  is  the  behavior  of  the 
energy  derivative  terms  which  are  causing  the  problem.  A  look  into 
the  equation  of  state  tables  for  the  material  of  this  zbne  showed  a 
discontinuity  in  these  values  near  this  temperature  and  density.  As 
it  turned  out,  a  small  increase  in  temperature  resulted  in  a  huge 
change  in  the  derivatives  due  primarily  to  the  linear  interpolation 
scheme.  Subsequently,  the  temperature  dropped  and  again  the  corres¬ 
ponding  derivatives  changed  drastically.  Thus,  the  temperature  (and 
corresponding  quantities  affected  by  it)  oscillated  back  and  forth 
about  an  entry  in  the  equation  of  state  tables,  either  side  of.  Which 
gave  widely  varying  derivative  values  based  upon  linear  interpolation. 
Decreasing  and  from  0.5  to  0.1  dampened  out  the  oscillations  and 
the  calculations  were  able  to  proceed.  Alternatives  Would  be  to  raise 
the  convergence  limit,  reduce  the  time  step  Size,  or  any  combination 
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Figure  28 
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CONCLUSION  AND  AREAS  FOR  FUTURE  RESEARCH 


As  stated  in  the  introduction,  the  development  of  this  system  hae 
extended  over  the  past  several  years.  During  this  period  a  variety  of 
techniques  have  been  tried  within  .  number  of  different 

aoletieae  on  four  different  computers.  Ihe  end  produce  co„slst8  of  a 
hfShly  PdUshed  machine  language  program  which  consume.  huge  amounts 

tlTO  the  PDP-10  *  ™  to  Che  BnaU  amounts  of 

tiae  consumed  ,  the  relative!*  inefficient  and  unfiniehed  FORTRAN 
version  r>n  the  CDC-7600. 

The  PDP-10  veraion  was  orlgi„ally  coded  in  FORTRAN  only  to  find 
out  that  the  FORTRAN  system  on  that  computer  was  full  of  gross  error. 

e  ficiencies.  To  obtain  significant  results  on  the  PDP-10,  the 

calculations  for  4-3^  and  av  « 

3  and  AV  necessitated  using  double  precision. 

After  quickly  deteraining  that  the  DRC-.upplied  double  precision 
routines  were  overly  restrictive  on  the  range  of  operand,  with  which 
they  could  produce  significant  results  (and  contained  ainor  error.) 
the  aethod  was  reprograaaed  in  MACRO-10,  the  aachine  language  for  the 
PDP-10.  As  alluded  to  above,  this  also  included,. -couplets  package 
of  library  routines  for  perforulng  double  precision  arithaetic 
(addition,  subtraction,  multiplication,  and  division)  in  addition  to 
sons  elementary  mathematical  functions  (square  root,  natural  log  and 
exponential).  To  say  the  least,  it  i,  .ad  to  think  that  a  computer 
Of  this  sise  has  no  double  word  arithaetic  capability  and  that  none 
of  the  modern  computer,  heve  the  hardware  capability  to  perform  the 
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elementary  mathematical  functions.  It  m3y  be  of  interest  at  this 
point  to  mention  that  the  double  precision  square  root  routine  within 
the  PDP-10  vers  ten  of  HYRAD1  consumes  over  60X  of  the  total  calcula- 
tional  time. 

In  addition  to  the  fact  that  the  MACRO-10  version  works  and  the 

FORTRAN  version  didn't*  the  resultant  speed  and  storage  improvement 

of  the  MACRO-10  version  over  the  FORTRAN  version  on  the  PDP-10  was 

dramatic.  The  differences  between  the  PDP-10,  CDC-6600  and  CDC-7600 

ate  even  more  outstanding.  A  particularly  simple  problem  which  ran 

for  64  minutes  on  the  PDP-10  in  29k  words  of  core,  used  but  6  minutes 

on  the  CDC-6600  and  a  mere  minute  on  the  CDC-7600  in  32k  words  of 

core.  Remember  that  this  comparison  is  between  a  very  efficient 

"hand  crafted"  program  on  the  PDP-10  with  a  FORTRAN  version  on  the 

CDC  machines  that  is  well  known  for  its  inefficient  use  of  machine 

resources.  It  is  worth  noting  that  the  same  problem  consumed  .4 

3 

minutes  using  the  non-iterative  F  Los  Alamos  program  and  that  HYRAD1 
used  an  average  of  6.1  Iterations  per  time  step,  both  taking  nearly 
2900  time  steps.  It  is  not  unusual  for  a  problem  to  run  for  several 
hours  on  the  CDC-7600  computer. 

It  is  clear  that  increased  hardware  functions  are  necessary  for 
numerical  calculations  of  this  type.  Not  fast  simple  parallel  functions 
as  exhibited  by  the  ILLIAC  IV  computer,  but  rather  independently  func-  . 
tloning  units  with  increased  capability  (l.e.  square  root,  exponential, 
etc.).  Raw  computing  power  As  not  enough  though.  A  better  man-machine 
interface  needs  to  be  developed.  The  user  needs  to  retain  control  over 
the  calculations  rather  than  submitting  to  the  crude  demands  of  the 


machine . 
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A*  mentioned,  this  prog.1*:?  was  developed  under  a  simple  inter¬ 
active  (and  sometimes  graphical)  system.  The  user  was  thereby  able 
to  alter  parameters,  intioduce  data  changes,  direct  program  control, 
and  even  make  crude  prograr  patches.  This  was  done  to  get,  first  of 
all,  answers  which  may  otherwise  be  very  difficult  to  obtain  under  the 
traditional  "batch"  mode  of  processing.  Secondly,  through  the  inter¬ 
active  mode  of  operation,  the  user  was  able  to  get  better  answers  in 
the  sense  that  more  resolution  could  be  obtained  at  those  times  when 
it  was  required.  Finally  it  was  possible  to  get  an  immediate  and 
in  depth  realization  of  the  solution  in  moving  picture  form. 

Slowly,  but  most  certainly,  tools  are  being  created  to  provide 
such  a  man-machine  interface.  However,  such  tools  have  traditionally 
been  improved  interface  equipment  such  as  display  consoles  or  they 
have  been  huge  and  rather  obtrusive  software  systems  which  are  programs 
to  make  up  for  the  lack  of  sufficient  and  sophisticated  hardware. 

Very  large  compilers  are  continually  being  developed  at  great  expem* 
to  provide  a  more  natural  language  interface  between  man  and  machine. 
However,  all  of  them  fall  short  on  retaining  enough  information  to 
relate  bSdk  to  the  original  language  and  to  retain  an  overall  view  of 
what  is  being  done.  Thus,  grand  interpretive  systems  have  been  and 
are  being  developed  which  have  some  of  these  characteristics  at  great 
expense  in  both  time  and  3pace. 

Clearly  the  developr  uit  of  au^h  software  syatems  over  the  past 
several  years  is  causing  questions  tc  be  asked.  Questions  such  as: 

Why  not  build  machines  to  work  at  the  u.^er  level?  Why  have  compilers? 


Thty  generally  discard  most  all  of  the  useful  information  required  os 
input  anyway.  Hardware  is  getting  so  fast  and  cheap,  why  not  build 
computara  that  work  in  tha  infix  mode  directly  (essentially)  on  the 
user's  source  program  which  can  then  be  easily  changed  and  debugged 
dynamically?  Thus,  why  not  build  computers  which  add,  multiply,  ate. 
operands  according  to  their  definition  at  the  point  of  execution  be 
it  integer,  real,  array,  procedure  or  whatever? 

There  era  so  many  situations  which  arise  at  run  time  that  Just 
can't  be  taken  care  of  by  the  compiler,  or  the  prograoner  ahead  of 
time.  This  is  true  of  problems  in  general  despite  it.a  best  planning. 
Therefore,  why  not  defur  the  find  decision  making  until  such  a  time 
^  context?  This  technique  is  already  in  popular  use  for 

detarmining  tha  final  operand  address.  It  saeaa  to  be  a  rather 
natural  extension  to  go  one  step  further  *»rd  eay  that  the  add  instruc¬ 
tion,  for  example,  la  defined  by  the  type  and  kind  of  operands  upon 
which  it  is  to  oparate.  This  than  makes  "add"  a  primitive  which  has 
no  complete  or  definitive  meaning  by  iteelf.  Thue,  the  actual  instruc¬ 
tion  set  for  a  computer  becomee  very  small  and  slapla.  There  are  no 
longer  2-5  different  multiply  instruct! ,ns  plu-  a  host  of  multiply 
subroutines  to  handle  things  like  the  multiplication  of  arrays  or 
complex  numbers.  Thara  is  just  one  multiply,  and  it  is  micro-programed, 
i*  you  wish,  by  its  operands. 

In  conclusion,  tha  model  developed  heroin  and  the  subsequent  cal¬ 
culations  indicate  that  the  numerical  methods  era  reasonable  and  are 
ca  ’able  of  giving  valid  results.  With  tha  advent  of  batter  computing 
tools,  the  solution  techniques  appear  even  more  palatable.  The  develop- 
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Met  of  the  straight  forward  radiation  diffusion  schese  was  an  ssrly 
(5  years  ago)  revolt  against  the  firm  order  expansion  techniques  as 
Ufesd  by  [5].  The  results  seem  to  bear  ou  the  fact  that  toese  equations 
are  at  leest  as  good  as  thos^  used  In  tho  pact.  Other  recent  enalyais 
and  work  done,  notably  by  Carton  Wendroff  [15],  soem  to  Indicate  thet 
this  approech  has  a  gre>t  deal  of  promise. 

The  Iterative  procedure  was  Included  for  two  purposes  First,  It 
provides  e  corrective  mechanism  which  tends  v,  give  better  than  flrat 
order  resolution  since  It  Includes  sons  foresight  ««  veil  as  hindsight. 
Secondly,  it  providts  an  easy  way  of  Incorporating  non-linear  effects. 

It  is  actually  Instructive  to  watch  the  iterative  procedure  sneak  up 
to  the  answer  vi*  the  dl*,>lav  console.  The  anomalies  ere  eepeclelly 
interesting  as  pointed  out  in  section  7.6.  Such  Iterative  procedures 
have  >een  successfully  used  for  a  long  tine  tu  thermal  and  neutron 
diff islon  programs 

Several  othir  foatures  of  the  model  are  also  worth  Mntloning  at 
this  time.  Probably  next  in  order  of  Importance  Is  the  opacity 
averaging  scheme  at  the  Interface  as  given  In  section  5.2  even  though 
its  worth  has  not  yet  be«n  fully  determined.  It  appears  to  be  the 
only  scheme  in  use  that  has  e  sound  mathematical  and  physical  basis 
and  yat  still  works  [6],  This  has  bee'.i  and  continues  to  be  a  prlM 
araa  for  future  work.  Coming  nest  In  a  close,  yet,  second  place  la 
tha  tons  centering  problem  as  exhibited  by  the  point  source  end  fruttnas 
calculations  in  section  7.3.  The  kinetic  energy  as  approximated  by 
(2.66)  or  the  total  energy  conservation  Is  Incorrect  for  these  problems. 
This  discrepancy  seems  to  be  conntcted  with  the  tone  centering  problem. 


107 


It  It  clear  that  tort  work  nee  a  to  be  dona  In  enla  araa  [16].  Evan  the 
3 

F  solution  of  tha  source  problea  la  In  suspect  In  light  of  the  ahock 
tubt  probla: .  it  seems  as  though  the  cslculstlonol  rtaulta  should  lsg 
bthlnd  sad  achiave  an  asymptotic  solution.  Problem  with  the  elngular- 
Itlaa  at  tha  canter  plus  the  fact  that  the  source  energy  la  actually 
Introduced  at  t*is  caucr  of  the  first  rone  rather  than  actually  at  tha 
canter  of  the  sphere  tend  to  suggest  a  much  harder  look  at  r'ia  p rob lea 
end  tha  Idealised  solutions  given  In  [12],  fi3j  and  [14]. 

Soaa  of  the  other  new  and  successful  features  Include  the  Hydro¬ 
dynamic  Tins  Step  Restriction  procedure  vhlch  guarantees  that  not  only 
vUl  aonaa  not  cross,  but  that  they  will  aot  expand  or  coapraaa  too 
quickly.  The  Implementation  of  the  extensive  equation  of  state  and 
opacity  tabla  lookup  procedural  are  not  detailed  herein,  but  era 
extensive,  efficient  and  vary  fast.  Tha  tachnlqua  of  Introducing 
OMrqy  Into  thn  problaa  by  aeons  of  temperature  profile  over  several 
sonan  as  given  la  aectlon  6.2  la  new  In  that  It  calculates  tha  amount 
of  onorgy  required  to  bring  tha  specified  tones  to  that  predetermine 

temperature.  Thn  Iterative  procedure  plays  a  useful  corrective  role 
hare  too. 

Tha  one  leportant  feature  not  Incorporated  within  tha  modal  la  a 
dynamic  reaonlng  procedure  iver  which  the  user  has  final  control,  yet 
ie  automatic  un^ar  certain  user  specified  conditions.  Thn  F3  program 
employe  a  procedjra  which  will  im  lnlt.'rlly  calculate  those  exterior 
•once  three  or  nora  aonaa  away  from  tha  #<tlve  senes.  This  technique 
help*  to  reduce  the  cslculstloe  tine  by  not  Including  thoee  tones 
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which  art  not  yat  active  (i.%.  Ooee  tenet  which  have  not  yet  bean 
affected  by  radiation  or  hydrodyaislc*).  ThU  »<-  <«*•  work*  only  f r* 
tha  so  called  ‘'Inner  boundary*'  (low  tone  nuebere)  outward  and  ouet 
hav*  an  ''aablent'*  boundary  condition  at  the  other  and  If  the  'outer 
boundary  condition  la  active.  If  there  le  a  gradient  near  the  "outer** 
boundary,  or  If  there  are  non-*ert>  eource  terra  near  the  "outer"  bouwd 
ary,  then  the  p-ogran  «u*t  proceed  to  do  the  calculation  for  all  aonaa 
rtfardlaaa  of  what  happen*  at  the  "Inner*  boundary  or  interior  to  the 
problen.  It  le  deelreble  to  general lie  the  technique  with  re toning 
into  wore  to nee  when  needed,  end  llkevlee,  reeonhlnlng  aonee  whan  no 

longer  needed. 
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